Mediation using maximum likelihood estimation with robust standard errors

[MP]

Hi all, I am trying to reproduce a SEM mediation analysis from a published article on my own data. The authors reported that parameter estimates were calculated using "maximum likelihood with robust standard errors" and that the "missing data were handled using full information maximum likelihood procedures". In an attempt to reproduce this, I entered this input:

modelFear <- 'CLg ~ c* HRR
             Fear ~ a* HRR
             CLg ~ b* Fear
             ab := a*b
             total := c + (a*b)'
fitFear <- sem(modelFear, data = cabsfull, missing = "fiml", se = "robust.sem")

When I run this, I get the following message:

Warning message:
In lav_options_set(opt) :
  lavaan WARNING: missing will be set to “listwise” for se = “robust.sem”

Does mean I cannot use full information ML to deal with missing values if I want robust standard errors?

If I omit the se argument, and input just the following:

fitFear <- sem(modelFear, data = cabsfull, missing = "fiml")

I get the following message:

Warning message:
In lav_data_full(data = data, group = group, cluster = cluster,  :
  lavaan WARNING: 5 cases were deleted due to missing values in 
		  exogenous variable(s), while fixed.x = TRUE.

So if I understood correctly, it used listwise deletion even though I asked for full information ML. Is it possible to reproduce what these authors have done - i.e. calculate the parameter estimates using the ML with robust standard errors and use FIML for missing data? If yes, what am I doing wrong?
I am new to R and SEM so any advice is very welcome. Thanks!

[Terrence Jorgensen]

You need Huber-White SEs (also a robust Yuan-Bentler test statistics).  You can set fixed.x=FALSE, or add ".x" to "fiml".  Try 

fitFear <- sem(modelFear, data = cabsfull, missing = "fiml.x", estimator = "MLR")

Terrence D. Jorgensen

Assistant Professor, Methods and Statistics

Research Institute for Child Development and Education, the University of Amsterdam

http://www.uva.nl/profile/t.d.jorgensen

[MP]

Hi Terrence, thank you very much for your reply! The solution you provided worked perfectly!

I was just wondering if I may ask you another question. As I mentioned before, I am trying to reproduce a mediation analysis that another group published some time ago. They mentioned  that "parameter estimates were calculated using maximum likelihood with robust standard errors" to account for non-normality in the data (my data also has the same issues) and "to test the significance of the indirect effects, bias‐corrected confidence intervals (CIs) for the indirect effects were generated using 1,000 bootstrap samples".  Therefore, If I use the function:

fitFear <- sem(modelFear, data = cabsfull, missing = "fiml.x", estimator = "MLR")




That will reproduce everything they did except bootstrapping CIs (bias-corrected) for the indirect effect. I was wondering if there is a possibility to add another function after my fitFear function that will just bootstrap the confidence intervals, or do I have to change my function entirely (and maybe try something like:





fitfear<- sem(ModelFear, data=cabsfull, missing = "fiml.x", estimator = "ML", test = "Yuan.Bentler", se = "bootstrap", bootstrap = 1000)

The overall idea is to provide the best solution that will account for the non-normality of the data. My sample size is N = 487. Thank you!

[PD]

You can only bootstrap if you have the raw data otherwise use the robust standard errors you have calculated to see if you get to the same conclusion that the authors did. If you want to bootstrap for your study you can use se = "bootstrap" and bootstrap = 1000 (although you can use 5000 if you want, it shouldn't make too much difference).

[Terrence Jorgensen]

I was wondering if there is a possibility to add another function after my fitFear function that will just bootstrap the confidence intervals

 

Use the parameterEstimates() function, which has an argument boot.ci.type="bca.simple" to request BC bootstrap CIs

[MP]

Thanks Terrence! Just a final question, how do I enter a number of bootstrap draws in the parameterEstimates() function? I know that for the sem() function it is bootstrap = , but lavaan does not recognise that argument if I pass it to the parameterEstimates() function. If I just enter:

parameterEstimates(fitFear, boot.ci.type="bca.simple")

It works fine, but it doesn't say on how many bootstrap samples it is based on.

Thanks for all the help!

[Terrence Jorgensen]

parameterEstimates(fitFear, boot.ci.type="bca.simple")

It works fine, but it doesn't say on how many bootstrap samples it is based on.

Because you already specified that when you fit the model fitFear

[MP]

Apologies, but I am a bit confused now, because my fitFear model does not include a bootstrapping argument. It is as you recommended earlier:

fitFear <- sem(modelfear, data = cabsfull, estimator = "MLR", missing = "fiml.x")                                                                        parameterEstimates(fitFear, boot.ci.type="bca.simple")

(The first line should give me an MLR estimator with robust standard errors and FIML for missing data, and the second should just provide me with BCa bootstrapped CI for indirect effect)            

Also, when I run both fitFear and parameterEstimates functions I get the results immediately, whereas when I usually bootstrap things in Lavaan it takes time as it's requires computational power, so I'm not sure if I am doing it properly?

Many thanks for your patience!

[Terrence Jorgensen]

Apologies, but I am a bit confused now, because my fitFear model does not include a bootstrapping argument.

Then add it.

It is as you recommended earlier:

You mean in response to your first post, which did not mention bootstrapping at all?

[MP]

Thanks! Just to confirm, it should be like this then?

fitFear <- sem(modelfear, data = cabsfull, estimator = "MLR", missing = "fiml.x", bootstrap = 5000)                                            parameterEstimates(fitFear, boot.ci.type="bca.simple")

[Terrence Jorgensen]

On Monday, April 29, 2019 at 3:11:33 PM UTC+2, MP wrote:

Thanks! Just to confirm, it should be like this then?

fitFear <- sem(modelfear, data = cabsfull, estimator = "MLR", missing = "fiml.x", bootstrap = 5000)                                            parameterEstimates(fitFear, boot.ci.type="bca.simple")

Sure, but I don't see why you would insist on MLR for robust SEs when you are already bootstrapping.  And if your N is large, there isn't any point to bootstrapping, just rely on the delta-method (robust) SEs provided by default for your indirect effects.

[Nikola Ćirović]

Hello!

I also have a question regarding this. Namely, I followed the example above and specified the following command:

fit.<-sem(model,  data = data ,orthogonal=TRUE, estimator='ML', missing = "fiml.x", estimator = "ML", test = "Yuan.Bentler", se = "bootstrap", bootstrap = 500)

However, after this, I was confused to find that the scaled fit indices (CFI/TLI) have risen substantially (robust ones did not though) in comparison to what is obtained with the following code:

fit<-sem(model,  data = data,std.lv=TRUE,orthogonal=TRUE, estimator='MLR', missing='fiml')

I am not sure if the command also bootstrapped the chi-square statistic or only the SE and CI (the CI is important to me as this is a latent mediation model). The scaled chi-squared after first (bootstrapped) and second code differed only slightly (YB scaling factor differed in third decimal) and it confuses me that CFI is much more favourable after bootstrap. 

I understand the recommendation to simply use the delta method, however, I am writing in hope to understand what the command I have written actually performed. 

Sorry if I am asking (missing) something basic. 

Thanks in advance.

[Terrence Jorgensen]

On Sunday, May 17, 2020 at 12:36:38 AM UTC+2, Nikola Ćirović wrote:

Hello!

I also have a question regarding this. Namely, I followed the example above and specified the following command:

fit.<-sem(model,  data = data ,orthogonal=TRUE, estimator='ML', missing = "fiml.x", estimator = "ML", test = "Yuan.Bentler", se = "bootstrap", bootstrap = 500)

However, after this, I was confused to find that the scaled fit indices (CFI/TLI) have risen substantially (robust ones did not though) in comparison to what is obtained with the following code:

fit<-sem(model,  data = data,std.lv=TRUE,orthogonal=TRUE, estimator='MLR', missing='fiml')

I am not sure if the command also bootstrapped the chi-square statistic or only the SE and CI (the CI is important to me as this is a latent mediation model). The scaled chi-squared after first (bootstrapped) and second code differed only slightly (YB scaling factor differed in third decimal) and it confuses me that CFI is much more favourable after bootstrap. 

Setting estimator = "MLR" sets test="yuan.bentler.mplus", not test="yuan.bentler".   See ?lavOptions description.  The Mplus variant is only asymptotically equivalent to what Yuan & Bentler described.

If the result didn't change much for the hypothesized model, then CFI probably changed because the baseline model's chi-squared changed more substantially.

FYI, setting missing = "fiml" vs. "fiml.x" can yield different behavior too, if you have missing data on exogenous observed variables, "fiml" will still listwise-delete observations with missing data on the exogenous covariates, but "fiml.x" retains them (so the model will be fitted to a larger sample).

[Nikola Ćirović]

Dear professor Jorgensen, 

Thank you so much for your quick reply. 

I understand that there may be more than one possible reasons for this (I don't have observed exogenous variables as far as I understand because my model is a measurement model with latent paths for which I want to bootstrap SE and get CI so I am trying to wrap my mind around this other reason you singled out as a possibility - the MLR). 

However, given the task at hand - and that is to obtain bootstrapped CI for the latent path coefficients, i.e. their indirect effects - I want to be safe and only take the se="bootstrap" option (only using default ML, even though I need MLR, however, once I obtain MLR fit indices it remains only to find the CI of the paths). It takes quite some time to do this (setting it to 500 boot takes an hour) and I judge that it really performs boot on all the SEs. However, I find it confusing that when I chose parameterEstimates function with "boot.ci.type = "bca.simple" argument (as I am trying to follow Heyes on this one and to use bias-corrected percentile bootstrap) - this gets finished in a second. So, I wonder how may I described what happened with these two commands. Can I use the same terminology that Heyes uses in his book on the process analysis and can I say that I obtained bias-corrected bootstrapped CI for the latent path coefficients - indirect and total effects?

[Terrence Jorgensen]

to obtain bootstrapped CI for the latent path coefficients, i.e. their indirect effects - I want to be safe and only take the se="bootstrap" option (only using default ML, even though I need MLR, however, once I obtain MLR fit indices it remains only to find the CI of the paths). It takes quite some time to do this (setting it to 500 boot takes an hour) and I judge that it really performs boot on all the SEs.

Yes, the entire model is fitted to each bootstrap sample.

 

However, I find it confusing that when I chose parameterEstimates function with "boot.ci.type = "bca.simple" argument (as I am trying to follow Heyes on this one and to use bias-corrected percentile bootstrap) - this gets finished in a second. 

Because the bootstrapping already happened when you called sem(), so parameterEstimates() only needs to do the calculations, not the bootstrapping. 

FYI, since you are using MLR anyway, it would be faster to use the Monte Carlo approach.  That would also take seconds rather than hours for 2000+ bootstraps (which is the condition under which bootstrapping works well, not 500), and performs just as well as bootstrapping (see References and Examples on the ?semTools::monteCarloMed help page).

[Patrick (Malone Quantitative)]

Chiming in that the principal reason to use bootstrapped CIs for indirect effects is to allow the CIs to be asymmetric, because the sampling distribution of a product of coefficients is asymmetric. Using bootstrap just to generate the SEs, while it has other merits in some situations, does not achieve this. Hence using the various approaches to boostrapped CIs (percentile, bias-corrected, accelerated).

--
You received this message because you are subscribed to the Google Groups "lavaan" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lavaan+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/lavaan/4a962340-269b-434a-afce-6f7ee49a9a12%40googlegroups.com.

--

Patrick S. Malone, Ph.D., Malone Quantitative
NEW Service Models: http://malonequantitative.com

He/Him/His

[Nikola Ćirović]

Dear Mr. Malone - thank you for your additional clarification. 

Dear Mr. Jorgensen, thank you for feedback. 

If I understood you correctly the procedure I used (the code for bias-corrected CI) is well written but the very procedure is insufficient in my model given the circumstances. 

I decided to pursue your advice to use Montecarlo method but it is somewhat more difficult to use that function (however I am eager to explore semTools' capabilities). Could you please help me by checking whether I used the function correctly - it is in the bottom of my code for one of the effects (I copied the syntax as it also confuses me how to correctly mark the parameters - so I use the same letters as I used in the syntax).

I used some examples for this google group. 

model.temp.rum<-'   

bis=~ p114+p125+p129+p130+p133+p136+p138

bas=~p116+p118+p127+p131+p137+p140

fight=~p115+p121+p123+p124+p132+p139

freez=~p117+p119+p126+p128+p134

flight=~p113+p120+p122+p135+p141

bis~~bas

bis~~fight

bis~~freez

bis~~flight

bas~~fight

bas~~freez

bas~~flight

fight~~freez

fight~~flight

freez~~flight

SAA=~ p142+p143+p144+p145+p146+p147+p148+p149+p150+p151+p152+p153+p154+p155+p156+p157+p158+p159+p160+p161+p162+p163+p164+p165+p166  

SA1=~   p143+p144+p145+p146

SA2=~  p144+p154+p155

SA3=~ p156+p157+p158+p159+p160

SA4=~p161+p162+p163+p164+p165+p166 

XY =~ p167+p168+p169+p170+p171+p172+p173+p174+p175+p176+p177+p178+p179+p180+p181+p183+p184+p185+p186

bbp =~p167+p168+p169+p170

vdp =~ p172+p173+p174

fft    =~p175+p177+p178+p179

p171 ~~ p185

p185 ~~ p186

XY~ b*bis+fight+flight+freez+a*bas

SAA~ c*XY+e*bis+fight+flight+freez+d*bas

direktBAS:= d

direktBIS:= e 

indirektBAS:=a*c

indirektBIS:=b*c

totalBAS:=d+(a*c)

totalBIS:=e+(b*c)

'

fit.model.temp.rum <-sem(model.temp.rum,  data = data ,std.lv=TRUE,orthogonal=TRUE, estimator='MLR', missing = "fiml")

summary(fit.model.temp.rum,

        fit.measures = TRUE,standardize = TRUE,

        rsquare = TRUE)

# only for indirect effect BAS-XY-SAA

med <- 'a*c'

myParams <- c("a","c")

myCoefs <- coef(fit.model.temp.rum)[myParams]     

myACM <- vcov(fit.model.temp.rum)[myParams, myParams]         #This one I did not particularly understand 

monteCarloMed(med, myCoefs, ACM = myACM, rep = 2000)

This produced a more narrow interval than my bias-corrected CI with 500 reps. 

Thanks.

[Terrence Jorgensen]

the principal reason to use bootstrapped CIs for indirect effects is to allow the CIs to be asymmetric

Likewise, Monte Carlo CIs make no assumption about the form of the distribution of the product (or any function of model parameters), so can be asymmetric. See the references in the ?monteCarloMed help page.

 

med <- 'a*c'

myParams <- c("a","c")

myCoefs <- coef(fit.model.temp.rum)[myParams]     

myACM <- vcov(fit.model.temp.rum)[myParams, myParams]         #This one I did not particularly understand 

monteCarloMed(med, myCoefs, ACM = myACM, rep = 2000)

You are combining 2 different examples from the help page.  Notice the comments say the second example is unnecessary.  Mimic the first lavaan example:

monteCarloMed('a*c', object = fit.model.temp.rum)

[Nikola Ćirović]

Dear Mr. Jorgensen, 

Thank you! This works well. So the single line of command that you wrote in your last post is enough? (I was confused with the command that specifies the "Asymptotic covariance matrix" )

I just have one more question (again maybe I'm missing something obvious) - can I use this command for total effects too? (the help page mentions only indirect effect - products - so it left me confused about the total effects).

Like this: 

monteCarloMed('d+(a*c)', object = fit.model.temp.rum.SA.basicBEZantiZAEFEbezBSre)

[Terrence Jorgensen]

So the single line of command that you wrote in your last post is enough? (I was confused with the command that specifies the "Asymptotic covariance matrix" )

 

Yes.  When you pass a lavaan model, it will detect that and automatically extract the coefficients and their ACOV from the lavaan object.

 

can I use this command for total effects too? 

Like this: 

monteCarloMed('d+(a*c)', object = fit.model.temp.rum.SA.basicBEZantiZAEFEbezBSre)

Yes!  Honestly, it can be used to test any null hypothesis about a transformation of parameters.  The Monte Carlo method can be used as a semi-parametric alternative to the generic, normal-theory-based "delta method" (e.g., ?deltaMethod in the car package).  Normal theory is assumed for individual estimated parameters, but (unlike the delta method) no parametric assumptions are made about any functions of parameters.

[Nikola Ćirović]

Dear Prof Jorgensen, thank you for your kind and helpfully elaborated comments (as always!).