bootstrapLavaan not working with standardizeSoultion

[Jes Coyle]

Upon investigating further my previous post (https://groups.google.com/forum/#!topic/lavaan/7w4m2EeVM38) I think the easiest way to calculate bootstrapped confidence intervals on the standardized coefficients is to use the bootstrapLavaan() function coupled with the standardizedSolution() function. However, when I try this with the example data it doesn't work. See code below:

HS.model <- '

visual =~ x1 + x2 + x3

textual =~ x4 + x5 + x6

speed =~ x7 + x8 + x9 '

fit <- cfa(HS.model, data=HolzingerSwineford1939, se='robust')

standardizedSolution(fit)

fit_std = bootstrapLavaan(fit, R=10, FUN=standardizedSolution)

I get the following error message:

Error in t.star[b, ] <- res[[b]] : 

  incorrect number of subscripts on matrix

Am I doing this correctly?

[Jes Coyle]

Fixed this by using the function

FUN=function(x) standardizedSolution(x)$std.est

[Jes Coyle]

Sorry, typo. It was:

FUN=function(x) standardizedSolution(x)$est.std

It would be great if there was a way to output vectors of results into an array rather than being limited to a single value to matrix.

e.g.

function(x) c(coef(x), standardizedSolution(x)est.std)

--
You received this message because you are subscribed to a topic in the Google Groups "lavaan" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/lavaan/pIXy7NUa3c0/unsubscribe?hl=en.
To unsubscribe from this group and all its topics, send an email to lavaan+un...@googlegroups.com.
To post to this group, send email to lav...@googlegroups.com.
Visit this group at http://groups.google.com/group/lavaan?hl=en.
For more options, visit https://groups.google.com/groups/opt_out.
 
 

[yrosseel]

On 04/18/2013 09:42 PM, Jes Coyle wrote:
> Fixed this by using the function
>
> FUN=function(x) standardizedSolution(x)$std.est

Ah, you figured it out yourself.

Yves.

[yrosseel]

On 04/18/2013 09:49 PM, Jes Coyle wrote:
> Sorry, typo. It was:
>
> FUN=function(x) standardizedSolution(x)$est.std
>
> It would be great if there was a way to output vectors of results into
> an array rather than being limited to a single value to matrix.
>
> e.g.
>
> function(x) c(coef(x), standardizedSolution(x)est.std)

This works for me?

example(cfa)

bootstrapLavaan(fit, R=10, FUN=function(x) c(coef(x),
standardizedSolution(x)$est.std), verbose=TRUE)
... bootstrap draw number: 1 OK -- niter = 31 fx =
0.160073283
... bootstrap draw number: 2 OK -- niter = 26 fx =
0.199542030
... bootstrap draw number: 3 OK -- niter = 27 fx =
0.203816957
... bootstrap draw number: 4 OK -- niter = 25 fx =
0.150926245
... bootstrap draw number: 5 OK -- niter = 24 fx =
0.161844035
... bootstrap draw number: 6 OK -- niter = 29 fx =
0.155509253
... bootstrap draw number: 7 OK -- niter = 25 fx =
0.158235483
... bootstrap draw number: 8 OK -- niter = 27 fx =
0.152239527
... bootstrap draw number: 9 OK -- niter = 24 fx =
0.182188231
... bootstrap draw number: 10 OK -- niter = 26 fx =
0.163480755
Number of successful bootstrap draws: 10
visual=~x2 visual=~x3 textual=~x5 textual=~x6 speed=~x8 speed=~x9
[1,] 0.6056449 0.7454157 1.132504 0.9044614 1.4404872 1.2113618
[2,] 0.3914412 0.7268797 1.061381 0.8902339 0.9953873 0.9830214
[3,] 0.5022559 0.5852260 1.031874 0.8998520 1.0354398 1.0573653
[4,] 0.5874887 0.6033834 1.192735 0.9793307 0.7952628 0.8401460
[5,] 0.5200911 0.6055870 1.080349 0.9195603 0.9759872 0.8652581
[6,] 0.4950064 0.5719694 1.234826 0.9564534 1.2462329 1.0514661
[7,] 0.3436045 0.4338989 1.139204 0.9490114 0.8893787 0.6654369
[8,] 0.4007381 0.6099811 1.172422 0.9276068 1.2445515 0.7837231
[9,] 0.6895950 0.7912108 1.175760 0.9278514 1.2217241 1.0398842
[10,] 0.5741523 0.6840139 1.077109 0.9680779 1.1477551 1.0300589
x1~~x1 x2~~x2 x3~~x3 x4~~x4 x5~~x5 x6~~x6 x7~~x7
[1,] 0.6436863 1.1834345 0.9985335 0.3679127 0.5026720 0.2441066 0.7944380
[2,] 0.3774312 1.2003884 0.9120006 0.3214846 0.4671205 0.4065354 0.7010751
[3,] 0.4462759 1.0496592 1.0213364 0.3200865 0.5288105 0.4156903 0.8381299
[4,] 0.3610213 0.8568051 0.8843278 0.3176423 0.4021452 0.2695766 0.6457401
[5,] 0.3900648 1.1901473 0.8922127 0.3213190 0.4602924 0.3716649 0.6373575
[6,] 0.5192102 1.0006174 0.8964864 0.3637627 0.3882180 0.2967474 1.1179113
[7,] 0.2080386 1.2674506 1.0681498 0.3344635 0.4599510 0.2508640 0.7013130
[8,] 0.4523860 1.1250614 0.8206319 0.3754805 0.3637526 0.3289645 0.7669970
[9,] 0.6463269 1.1347346 0.7272798 0.3608645 0.4966565 0.3762110 0.6768754
[10,] 0.4132028 1.0617631 0.7490325 0.3286215 0.4608382 0.3249921 0.6697479
x8~~x8 x9~~x9 visual~~visual textual~~textual speed~~speed
[1,] 0.5435101 0.6569349 0.6793346 0.9362303 0.2635697
[2,] 0.7160105 0.4592423 0.8967690 1.2142987 0.5904480
[3,] 0.4981861 0.5589465 0.9906092 1.1721806 0.4637479
[4,] 0.4344006 0.5569914 0.7817200 0.7246961 0.5411067
[5,] 0.5024131 0.5111072 0.9994318 0.9338323 0.5533503
[6,] 0.5090305 0.4628259 0.8228565 0.8943179 0.3086649
[7,] 0.4336595 0.5963139 1.2280367 0.9304618 0.6853922
[8,] 0.2829098 0.6376419 0.7692579 0.8230174 0.3992594
[9,] 0.4967879 0.5805021 0.6860882 1.1117231 0.4606628
[10,] 0.4283863 0.5957132 0.9440382 1.0892194 0.4907127
visual~~textual visual~~speed textual~~speed
[1,] 0.3245811 0.1907390 0.07376072 0.7165698 0.4170563
[2,] 0.5658651 0.3797039 0.34658075 0.8389218 0.3204880
[3,] 0.5826129 0.3090549 0.25891884 0.8303098 0.4385100
[4,] 0.2280253 0.1469083 0.15975986 0.8270879 0.4893711
[5,] 0.4471759 0.2819929 0.18506837 0.8481015 0.4302362
[6,] 0.4482657 0.2196899 0.16297120 0.7830239 0.4095215
[7,] 0.4946017 0.1818422 0.20330043 0.9247345 0.3203907
[8,] 0.4138314 0.1411709 0.10032015 0.7935306 0.3145469
[9,] 0.3574672 0.2171681 0.15705708 0.7175798 0.4725625
[10,] 0.4896993 0.3370024 0.26804895 0.8340005 0.4760932

[1,] 0.5237582 0.8472835 0.8395887 0.8708088 0.4991181 0.7082076 0.6087437
[2,] 0.5847237 0.8891966 0.8633929 0.8384614 0.6761448 0.6705641 0.7443465
[3,] 0.4993534 0.8862862 0.8380920 0.8339246 0.5968372 0.7067567 0.6937008
[4,] 0.4934291 0.8338225 0.8481685 0.8488447 0.6752182 0.6638163 0.6377931
[5,] 0.5396164 0.8625542 0.8384981 0.8245965 0.6817066 0.7155326 0.6690900
[6,] 0.4805566 0.8431245 0.8822687 0.8566378 0.4651533 0.6964225 0.6514623
[7,] 0.4218239 0.8576633 0.8509793 0.8772733 0.7030359 0.7453765 0.5807659
[8,] 0.5085188 0.8286781 0.8698823 0.8263254 0.5851006 0.8283226 0.5270361
[9,] 0.6093357 0.8688759 0.8693472 0.8472516 0.6363684 0.7619379 0.6795751
[10,] 0.6090516 0.8764839 0.8560208 0.8709250 0.6502771 0.7755226 0.6829229

[1,] 0.4865277 0.8260641 0.7256773 0.2821107 0.2950908 0.2416920 0.7508811
[2,] 0.2962103 0.8972875 0.6580983 0.2093294 0.2545527 0.2969824 0.5428282
[3,] 0.3105857 0.8077090 0.7506462 0.2144968 0.2976017 0.3045698 0.6437854
[4,] 0.3159256 0.7605160 0.7565277 0.3047400 0.2806102 0.2794628 0.5440804
[5,] 0.2807238 0.8148968 0.7088142 0.2560002 0.2969210 0.3200407 0.5352762
[6,] 0.3868736 0.8322921 0.7690653 0.2891410 0.2216020 0.2661717 0.7836324
[7,] 0.1448661 0.8973498 0.8220646 0.2644136 0.2758342 0.2303916 0.5057405
[8,] 0.3703092 0.9010602 0.7414087 0.3132926 0.2433047 0.3171863 0.6576573
[9,] 0.4850792 0.7766846 0.6287100 0.2450547 0.2442354 0.2821648 0.5950353
[10,] 0.3044432 0.7733353 0.6290561 0.2317760 0.2672284 0.2414896 0.5771397

[1,] 0.4984419 0.6294311 1 1 1 0.4069959 0.4507647 0.1484862
[2,] 0.5503438 0.4459483 1 1 1 0.5422629 0.5218118 0.4093086
[3,] 0.5004950 0.5187792 1 1 1 0.5406692 0.4559776 0.3511765
[4,] 0.5593479 0.5932200 1 1 1 0.3029559 0.2258808 0.2551222
[5,] 0.4880130 0.5523185 1 1 1 0.4628788 0.3791941 0.2574528
[6,] 0.5149958 0.5755968 1 1 1 0.5225496 0.4359176 0.3101853
[7,] 0.4444139 0.6627110 1 1 1 0.4627012 0.1982072 0.2545772
[8,] 0.3138817 0.7222329 1 1 1 0.5200953 0.2547308 0.1750073
[9,] 0.4194506 0.5381776 1 1 1 0.4093058 0.3862910 0.2194662
[10,] 0.3985646 0.5336164 1 1 1 0.4829220 0.4951361 0.3666425

[Jes Coyle]

Right- it works, it just cbinds the results rather than creates a 3-d array, which might be easier to work with as output. Not really a big deal, though.

--
You received this message because you are subscribed to a topic in the Google Groups "lavaan" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/lavaan/pIXy7NUa3c0/unsubscribe?hl=en.

To unsubscribe from this group and all its topics, send an email to lavaan+unsubscribe@googlegroups.com.

[JB]

Dear Jes and Yves,
Many Thanks for you very useful thread. I would like to make sure that I correctly extract the confidence intervals for standardized estimates, without and with boostrapping. It is very basic but I would be grateful if you could tell me if my syntax is correct. Building on your example:

fit <- cfa(HS.model, data=HolzingerSwineford1939, se='robust')

#confidence intervals of standardized solution (using the standard error given by standardizedSolution)
fit_std=standardizedSolution(fit)
results=data.frame(Names=paste0(parameterEstimates(fit)$lhs,parameterEstimates(fit)$op,parameterEstimates(fit)$rhs),
                   Estimate=fit_std$est.std,
                   LowCI=fit_std$est.std-1.96*fit_std$se,
                   HiCI=fit_std$est.std+1.96*fit_std$se)
results

#Bootstrapped estimates and Confidence Intervals, using quantile() to extract the 95%CI
fit_boot = bootstrapLavaan(fit, R=100, FUN=function(x) standardizedSolution(x)$est.std)
results_boot=data.frame(Names=paste0(parameterEstimates(fit)$lhs,parameterEstimates(fit)$op,parameterEstimates(fit)$rhs),
                     Estimate=colMeans(fit_boot),
                    LowCI=apply(fit_boot,2,FUN=function(x)quantile(x,0.025)),
                   HiCI=apply(fit_boot,2,FUN=function(x)quantile(x,0.975)))
results_boot                  
 
#compare results
data.frame(results,results_boot)

Two questions:
Regarding se='robust', I can understand that it is important for CI without boostrap. But when we use bootstrapping, I'm not sure why a robust standard error would make any difference as the confidence intervals are based on bootsrapping? 
Is it correct to use a mean to compute a boostrap estimate as I did with colMeans()
Many Thanks!
JB

[yrosseel]

On 01/14/2014 08:10 PM, JB wrote:
> Dear Jes and Yves,

> fit_std=standardizedSolution(fit)
> results=data.frame(Names=paste0(parameterEstimates(fit)$lhs,parameterEstimates(fit)$op,parameterEstimates(fit)$rhs),
> Estimate=fit_std$est.std,
> LowCI=fit_std$est.std-1.96*fit_std$se,
> HiCI=fit_std$est.std+1.96*fit_std$se)

Looks good. A similar (although with slightly different values) approach is

fit2 <- cfa(HS.model, data=HolzingerSwineford1939, se="robust",
std.lv=TRUE, std.ov=TRUE)
parameterEstimates(fit2)

giving you the lower/upper bounds directly.

> #Bootstrapped estimates and Confidence Intervals, using quantile() to
> extract the 95%CI
> fit_boot = bootstrapLavaan(fit, R=100, FUN=function(x)
> standardizedSolution(x)$est.std)
> results_boot=data.frame(Names=paste0(parameterEstimates(fit)$lhs,parameterEstimates(fit)$op,parameterEstimates(fit)$rhs),
> Estimate=colMeans(fit_boot),
>
> LowCI=apply(fit_boot,2,FUN=function(x)quantile(x,0.025)),
> HiCI=apply(fit_boot,2,FUN=function(x)quantile(x,0.975)))

Looks good.

> Two questions:
> Regarding se='robust', I can understand that it is important for CI
> without boostrap. But when we use bootstrapping, I'm not sure why a
> robust standard error would make any difference as the confidence
> intervals are based on bootsrapping?

It would not make any difference here. It depends on what you extract
from the bootstrapped model: here you extract only the 'est.std' values,
so it does not matter what the value of the se= argument was.

> Is it correct to use a mean to compute a boostrap estimate as I did with
> colMeans()

Yes!

Yves.

[JB]

Hello!
Following our previous discussion: in my models (quantitative genetic), the quantity of interest is the square of the std.all estimate rather than the estimate itself. I think I can compute a bootstrapped estimate by simply adapting the syntax above:

 fit_boot = bootstrapLavaan(fit, R=100, FUN=function(x) standardizedSolution(x)$est.std^2)
 
The ^2 at the end will do the trick, is that correct?

However, boostrapping is time consuming as one model takes almost one minute to run (lots of people, complex model). I'm not sure there are guidlines regarding the number of boostrapping but even 1000 would take a long time. So it would be good to obtain maximum likelihood confidence intervals of  est.std^2 . However, I'm unsure of how to do that... and as CI are more important than pvalue for this model, I would like to make sure I do that right...
Many thanks for your help!!!
JB
PS: I have convinced two more SEM users to switch to lavaan :)

[yrosseel]

On 02/07/2014 04:40 PM, JB wrote:
> Hello!
> Following our previous discussion: in my models (quantitative genetic),
> the quantity of interest is the square of the std.all estimate rather
> than the estimate itself. I think I can compute a bootstrapped estimate
> by simply adapting the syntax above:
>
> fit_boot = bootstrapLavaan(fit, R=100, FUN=function(x)
> standardizedSolution(x)$est.std^2)
>
> The ^2 at the end will do the trick, is that correct?

Yes.

> However, boostrapping is time consuming

True. If you have multiple cores on your machine, you can use the
parallel= options of the bootstrapLavaan() function to speed up things a
little. But lavaan can be slow, as it is currently written in 100% R.

> a long time. So it would be good to obtain maximum likelihood confidence
> intervals of est.std^2 . However, I'm unsure of how to do that...

Perhaps you could use the delta method, if you can work out the
derivative dy/dx, where x is the original parameter, and y is the
'transformed' parameter. The covariance of the original parameters can
be extracted by the vcov() function. But of course, this is just an
approximation.

Yves.

[JB]

 Dear Yves,
Thanks for your answer! I finally defined directly in my models all the quantities of interest I was interested in (e.g. res1:= label1*label2) and then got standard errors  and CI. If I am correct, the CI are based and standard errors and symmetric which is not always adequate for what I have to do.
In OpenMx, the mxCI function yields, maximum likelihood confidence intervals based on the article: "Michael C. Neale, Michael B. Miller, The Use of Likelihood-Based Confidence Intervals in Genetic Models, Behavior Genetics March 1997, Volume 27, Issue 2, pp 113-120". Does something like that exists in Lavaan?
Thanks again!!
JB

[yrosseel]

On 02/14/2014 12:42 PM, JB wrote:

> In OpenMx, the mxCI function yields, maximum likelihood confidence
> intervals

I actually got some research code from Carl F. Falk to do this. It
should have been included in lavaan a long time ago, but somehow, it is
still on my TODO list. Will bump it up.

Yves.

[JB]

That would be amazing!!! I was thinking of respecifying all my models in OpenMx with matrix specification to use mxCI, which would take me a considerable amount of time! Please keep us informed when you get round to it, and thank you so much for your work! 

[JB]

PS: that would be all the more interesting because on the few models that I run in both OpenMx and lavaan, lavaan was several time quicker. As the mxCI() usually take a lot of time so people sometimes give up trying to get the maximum likelihood CI. So lavaan would have an edge on that.