Specifying prior distributions

[Annemieke Brouwer]

Hi all,

I am trying to apply BSEM to estimate links between latent variables. To specify prior distributions, I have been looking at the example on covariance parameters on the website. It is stated that it is safest to stick with beta priors here, and I was wondering what the reason for this is. Could someone elaborate on this?

As prior input, I would like to have positive and negative expected relations between the latent variables. How would I formulate this using a beta distribution on the covariance parameters? The explanation on the website states blavaan automatically translates the prior to an equivalent distribution with support on (-1,1), but I've not been able to figure out what this means exactly. 

Thanks in advance!

Best,

Annemieke  

[Ed Merkle]

Annemieke,

There is a bit more information in Section 4.1 here:

https://arxiv.org/pdf/2008.07733.pdf

The short answer is that the prior distributions are on the correlations underlying the covariances. And we know that correlations are bounded by -1 and 1, which is where the beta distribution comes in.

But these prior distributions are impacted by the restriction that the correlation matrix stays positive definite, which means that the prior that you specify might end up being more informative than you think. There is some general detail about that issue here:

https://blog.mc-stan.org/2020/11/13/correlation-cats/

I would say that the priors on covariances and correlations are one of the trickiest parts of Bayesian SEM, and I am not completely happy with how blavaan currently handles them. I hope to improve this in the future.

Ed

--
You received this message because you are subscribed to the Google Groups "blavaan" group.
To unsubscribe from this group and stop receiving emails from it, send an email to blavaan+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/blavaan/2a1ad6c7-4c08-43eb-8e4c-4f6a388d13ccn%40googlegroups.com.

[Annemieke Brouwer]

Hi Ed,

Thanks a lot for your explanation and the links to further information on the topic.

If I understand correctly, a positive expected correlation could be formulated in the bsem model as follows:

academic_performance  ~~  prior("beta(1,1)") * autonomous_motivation

Is that correct? And would a beta(2,2) distribution mean that the probability of a correlation of 0.5 is highest? 

Lastly, would it be possible to express an expected negative correlation in the same way, through the beta distribution in prior()? What would the input be in the beta distribution?

academic_performance ~~ prior("beta(?,?)") * controlled_motivation

Thanks again for your time, 

Annemieke

Op woensdag 12 mei 2021 om 17:17:04 UTC+2 schreef Ed Merkle:

[Ed Merkle]

Annemieke,

The beta distributions for these correlations have support on (-1,1) instead of the usual (0,1), so your interpretations are a bit off:

beta(1,1) is uniform from -1 to 1.

beta(2,2) puts more density on correlations near 0, as opposed to extreme correlations.

To express the idea that positive correlations are more likely, you need something assymetric, like beta(5,3). And you could reverse those numbers if negative correlations are more likely.

To look at one of these distributions, you could do

x <- seq(.01, .99, .01)

plot(-1 + 2*x, dbeta(x, 5, 3), type="l")

Finally, related to my last email, these beta priors are usually more informative than the plot shows, because they must obey the constraint that the correlation matrix is positive definite.

Ed

[Mauricio Garnier-Villarreal]

Annemieke

As Ed mentions, with the correlations its harder to know how much information are your priors actually providing the model. In some cases, when I have added more strict priors, can become improper priors

Could you maybe explain in more detail what you are looking for to test? This way we might be able to guide a way to test it given the model setup

take care

[Annemieke Brouwer]

Hi Ed and Mauricio,

Thanks, the formulation of the beta distributions was a lot more clear after the explanation and the plot suggestion.

The model setup I’ve come up with now is stated below. I’m interested in the correlations between several latent variables, with answers to a survey as observed variables. For prior distributions, I’ve created a beta distribution with a correlation estimate as a mean and a variance of 0.01. The Bayes factor indicates that this model is better than the one without the priors, but it could be that its influence is more informative than wanted as Ed suggested. What do you think of this approach? And what would be the best way to compare these results to the results of the same model in Lavaan?

Thanks a lot!

Annemieke

library(blavaan)

library(readxl)

Data <- read_excel("~/3 Studie/BEP/Prime/R/Data/Data corrected.xlsx")

model<-'

#measurement part

SMC =~ SM1 + SM3 + SM4 + SM11 + SM12 + SM13 + SM14 + SM15 + SM16 + SM18 + SM19 + SM20 + SM21 + SM26

competence_satisfaction =~ Needs1 + Needs7 + Needs13 + Needs19

competence_frustration =~ Needs2 + Needs8 + Needs14 + Needs20

autonomy_satisfaction =~ Needs5 + Needs11 + Needs17 + Needs23

autonomy_frustration =~ Needs6 + Needs12 + Needs18 + Needs24

relatedness_satisfaction =~ Needs3 + Needs9 + Needs15 + Needs21

relatedness_frustration =~ Needs4 + Needs10 + Needs16 + Needs22

autonomous_motivation =~ Motivation2 + Motivation4 + Motivation7 + Motivation8 + Motivation11 + Motivation13 + Motivation15 + Motivation16

controlled_motivation =~ Motivation1 + Motivation3 + Motivation5 + Motivation6 + Motivation9 + Motivation10 + Motivation12 + Motivation14

academic_performance =~ Grade 

#structural part

relatedness_satisfaction ~ SMC

relatedness_frustration ~ SMC

competence_satisfaction ~ SMC

competence_frustration ~ SMC

autonomous_motivation ~ competence_satisfaction + competence_frustration + autonomy_satisfaction + autonomy_frustration + relatedness_satisfaction + relatedness_frustration 

controlled_motivation ~ competence_satisfaction + competence_frustration + autonomy_satisfaction + autonomy_frustration + relatedness_satisfaction + relatedness_frustration

academic_performance ~ autonomous_motivation + controlled_motivation

#priors

SMC ~~ prior("beta(48.69922,11.23828)") * competence_satisfaction #1

SMC ~~ prior("beta(48.69922,11.23828)") * relatedness_satisfaction #2

SMC ~~ prior("beta(33.47415,58.23585)") * competence_frustration #3

SMC ~~ prior("beta(33.47415,58.23585)") * relatedness_frustration #4

competence_satisfaction ~~ prior("beta(48.33765,10.97235)") * autonomous_motivation #5

autonomy_satisfaction ~~ prior("beta(48.33765,10.97235)") * autonomous_motivation #6

relatedness_satisfaction ~~ prior("beta(48.73498,11.265)") * autonomous_motivation #7

competence_satisfaction ~~ prior("beta(21.52807,57.04153)") * controlled_motivation #8

autonomy_satisfaction ~~ prior("beta(58.18198,33.769)") * controlled_motivation #9

relatedness_satisfaction ~~ prior("beta(21.20825,56.90685)") * controlled_motivation #10

competence_frustration ~~ prior("beta(30.83955,58.55045)") * autonomous_motivation #11

autonomy_frustration ~~ prior("beta(30.83955,58.55045)") * autonomous_motivation #12

relatedness_frustration ~~ prior("beta(30.83955,58.55045") * autonomous_motivation #13

competence_frustration ~~ prior("beta(56.2104,19.7496)") * controlled_motivation #14

autonomy_frustration ~~ prior("beta(56.2104,19.7496)") * controlled_motivation #15

relatedness_frustration ~~ prior("beta(56.2104,19.7496)") * controlled_motivation #16

academic_performance ~~ prior("beta(56.62967,38.94783)") * autonomous_motivation #17

academic_performance ~~ prior("beta(44.1,53.9)") * controlled_motivation #18

'

fit<-bsem(model, data=Data)

Op maandag 17 mei 2021 om 09:43:19 UTC+2 schreef mauga...@gmail.com:

[Mauricio Garnier-Villarreal]

Annemeike

If you will test a stronger hypothesis like this, would be recommended to

- estimate the model with different priors, so you can evaluate the effect of priors in general. You want the priors to add information to the model, but not to lead the results

- I dont recommend using the bayes factor for model comparison. It needs large posterior draws to be stable, and might present bias, that would need to be tested with simulations before its use in different setting https://arxiv.org/abs/2103.08744

- I would recommend you to use the LOO comparison. with blavCompare function, you get the LOO-IC for each model, and the comaprison between models with it. As long as you have a "large enough sample", this difference follows a normal distribution, so the general recommendation is to look at the ratio of LOO-dif/SE-diff, and then you can decide of a threshold to use to define is a model is "meaningfully better" than the other. The minimum recommneded is 2, which would be similar to the 95%, but its better to be more conservatiove with at least a ratio of 4 or highaer as the SE of the difference is optimistic https://arxiv.org/abs/2008.10296

- The LOO comparison is less sensitive to other model/priors characteristics as the Bayes Factor

- When using stricter priors, you whoucl plot some of these posterior distributions, to make sure that the priors are not giving strong bondaries to the estimates

- Another option, is the area inference approach. You estimate the model with less restrictive priors, and then you can set different ways to describe the posterior. For example, can present what proportion of the posterior is above a certain value like 0.3. If 90% of the posterior is above 0.3 for the correlation is strong evidence of the direction and magnitude of it. This way the focus is on describing areas of the posterior instead of setting string boundaries with priors

Hope this helps

[Annemieke Brouwer]

Mauricio,

Thanks a lot for your suggestions, I'm going to look into them. Very helpful!

I did run into some errors all of a sudden relating to the summary() function:

Error in getMethod("summary", "blavFitIndices") : 

  no method found for function 'summary' and signature blavFitIndices

I saw multiple threads in the group relating to the same issue, so I've tried the suggested options of putting library(lavaan) in front of library(blavaan); reinstalling the packages, I've even reinstalled RStudio, but the issues remain. I saw that Ed mentioned the extra function summary.blavaan(), but I couldn't find it (also not when using the github version of the package). Any ideas on what to try next to resolve the summary issues? 

As for code that causes the errors, it's these lines. 

1 AFIb <- ppmc(fitb, fit.measures = c("rmsea", "srmr", "cfi", "tli")) 

2 summary(AFIb)

3 blavFitIndices(fitb)

Line 3 gives an error, and line 2 gives the following: (instead of the posterior summary statistics) 

  Length    Class     Mode 

       1 blavPPMC       S4 

Annemieke

Op dinsdag 15 juni 2021 om 11:18:18 UTC+2 schreef mauga...@gmail.com:

[Ed Merkle]

Thanks for the report, I think I can fix this issue in the next day or so. I will let you know.

Ed

--

You received this message because you are subscribed to the Google Groups "blavaan" group.
To unsubscribe from this group and stop receiving emails from it, send an email to blavaan+u...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/blavaan/feb5be60-d495-4ab1-a1fc-1b7c4a52394fn%40googlegroups.com.

[Annemieke Brouwer]

Hi Ed,

Great, thanks a lot!

Annemieke

Op dinsdag 15 juni 2021 om 18:16:42 UTC+2 schreef Ed Merkle:

[Ed Merkle]

Annemieke,

I believe this issue is now fixed in the github version of blavaan. To install it, see the instructions near the bottom of this page:

https://github.com/ecmerkle/blavaan

Ed

[Annemieke Brouwer]

Ed,

Yes, it's working fine now, thanks a lot! 

Annemieke

Op woensdag 16 juni 2021 om 16:44:53 UTC+2 schreef Ed Merkle:

[Annemieke Brouwer]

Or at least, the summary function works fine now. For the BlavFitIndices I do get the following error:

> blavFitIndices(fitb)

Error in if (bopts$prisamp) { : argument is of length zero

Do you have an idea of what goes wrong here?

Best,

Annemieke

Op donderdag 17 juni 2021 om 12:20:22 UTC+2 schreef Annemieke Brouwer:

[Annemieke Brouwer]

Hi Ed,

If you have the time,  or someone else of course, could you help me with the following error that BlavFitIndices gives me currently? 

See my email below. If you need more information, let me know!

Best,

Annemieke

Op donderdag 17 juni 2021 om 13:15:24 UTC+2 schreef Annemieke Brouwer:

[Ed Merkle]

Sorry, I think I missed this message. I wonder whether your fitb comes from a previous version of blavaan, before you updated it last week. If so, you might try to re-estimate fitb to get it to work. Or, you might try this hack that I would not recommend as a general solution:

fitb@Options$prisamp <- FALSE

and then run blavFitIndices again.

Ed

--

You received this message because you are subscribed to the Google Groups "blavaan" group.
To unsubscribe from this group and stop receiving emails from it, send an email to blavaan+u...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/blavaan/3d5e5fd1-3b45-4421-a788-3497a83593a8n%40googlegroups.com.

[Annemieke Brouwer]

Hi Ed, 

Yes, after running the model again in the github version, everything works!

Thanks a lot. 

Op dinsdag 22 juni 2021 om 18:31:58 UTC+2 schreef Ed Merkle: