This section develops a Gibbs sampler to estimate SEM with reflective measurement indicators. ^{1 , 11 , 12} The Bayesian estimation is illustrated by considering a SEM that is equivalent to the mostly used model. A SEM is composed of a measurement equation (3) and a structural equation (4) ⁹ : y i = Λ ω i + ε i (3) η i = Π η i + Γ ξ i + δ i (4) where i ϵ 1 … n

It is assumed that measurement errors are uncorrelated with ω and δ , residuals are uncorrelated with ω and the variables are distributed as follows: ε i ∼ N 0 Ψ ε (5) δ i ∼ N 0 Ψ δ (6) ω i ∼ N 0 Σ ω (7)

∀ i ϵ 1 … n , where Ψ ε and Ψ δ are diagonal matrices. The covariance matrix of ω is derived based on the SEM: Σ ω = E ηη T E ξη T E ηξ T E ξξ T (8) Σ ω = Π 0 − 1 ΓΦΓ T + Ψ δ Π 0 − T Π 0 − 1 ΓΦ ΦΓ T Π 0 − T Φ (9) ηη T = Π 0 − 1 Γ ξ + Π 0 − 1 δ Π 0 − 1 Γ ξ + Π 0 − 1 δ T = Π 0 − 1 Γ ξξ T Γ T + δδ T Π 0 − T + Π 0 − 1 Γ ξδ T + δξ T Γ T Π 0 − T E ηη T = Π 0 − 1 ΓΦΓ T + Ψ δ Π 0 − T ηξ T = Π 0 − 1 Γ ξ + Π 0 − 1 δ ξ T E ηξ T = Π 0 − 1 Γ (10)

In order to enable Gibbs sampling from full conditional posterior distributions, natural conjugate prior distributions for the unknown parameters are considered. ²⁵ Let ψ εκ be the kth diagonal element of Ψ ε , ψ δι be the l th diagonal element of Ψ δ , Λ κ T be the kth row of Λ and M ι T be the lth row of M, ψ εk − 1 ∼ Gamma α 0 εk , β 0 εk (11) Λ k ψ εk − 1 ∼ N Λ 0 k ψ εk H 0 Λ k (12) ψ δi − 1 ∼ Gamma α 0 δi , β 0 δi (13) M i ψ δi − 1 ∼ N M 0 i ψ δi H 0 Mi (14) Φ ∼ IW v 0 V 0 (15) with κ ϵ 1 … p and ι ϵ 1 … q 1

The joint posterior of all unknown parameters is proportional to the likelihood times the prior, p Λ Ψ ε Ω M Ψ δ Φ Y ∝ p Y Λ Ψ ε Ω M Ψ δ Φ ∗ p Λ Λ ε Ω M Ψ δ Φ (16)

Given Y and Ω , Λ and Ψ ε are independent from Σ ω . Draws of Ω , can cause estimation of Λ and Ψ ε as a simple regression model. Thus, sampling from the posterior distribution of Λ and Ψ ε without reference to Σ ω . The same holds for inference with regard to M, Φ and Ψ δ , which are independent from Y given Ω .

The heteroscedastic error structure with different functional form of error variance under consideration are double logarithmic form, linear form, linear-inverse form and linear-absolute form as expressed in equation 17, 18, 19 and 20, respectively. σ 2 = ln σ 2 = λ o ∗ + λ i ∗ ln γ i ∗ + vi (17) σ 2 = | ε i ∗ ε i ∗ ' | = λ i ∗ + λ 2 ∗ γ i ∗ + νi 2 (18) σ 2 = | ε i ∗ ε i ∗ ' | = λ i ∗ + λ 2 ∗ γ i ∗ + νi 2 (19) σ 2 = | ε i ∗ ε i ∗ ' | = λ i ∗ + λ 2 ∗ γ i ∗ + νi 2 (20)

Each of the functional forms of heteroscedastic error structure will be incorporated into the modified model. The variance matrix for disturbance vector is given as ∑ = ε i ∗ ε j ∗ ' = σ λ 2 i ∗ i = j (21) Ω = σ 2 λ i ∗ 0 ⋯ 0 0 σ 2 λ 2 ∗ ⋯ 0 ⋮ 0 ⋯ 0 0 0 ⋯ σ 2 λ n ∗ (22)

The posterior density is the product of the likelihood and the prior distribution chosen ^{2 , 13} ( P λ ∗ h , Ω | y ∗ ) αp y ∗ λ ∗ h Ω p λ ∗ p h p Ω (23) p P λ ∗ h Ω y ∗ = h N 2 exp − h 2 γ ∗ − λ ∗ ' γ ∗ − λ ∗ γ ∗ × h N + v − k 2 exp − hv 2 s − 2 × p Ω × exp − 1 2 λ ∗ − λ 0 ∗ ′ V − 1 ¯ λ ∗ − λ 0 ∗ (24) P λ ∗ h Ω y ∗ = h N 2 | Ω | 1 2 exp − h 2 y ∗ − λ ∗ γ Ω − 1 y ∗ − λ ∗ γ × exp − v − 1 ¯ 2 λ ¯ ∗ − λ 0 ∗ V − 1 ¯ λ ¯ ∗ − λ 0 ∗ × n − 1 α + 1 exp − β h × Ω ∗ − 1 β 0 + k + 1 / 2 e − tr R 0 − 1 β 0 − 1 / 2 (25)

Since the full posterior distribution is intractable; a Markov chain Monte Carlo (MCMC) simulation method of Gibbs sampling is employed. ²⁵ This involves the use of marginal posterior distribution. λ = λ 0 ∗ = γ ′ Ω − 1 γ − 1 γ ' Ω − 1 γ = γ ∗ ' γ − 1 γ ∗ ' γ ∗ (26) S 2 = γ ∗ − γ ∗ λ 0 γ ∗ − γ ∗ λ 0 V ¯

Also Vs 2 ¯ + λ ∗ − λ ^ 0 ' γ i ∗ ' γ i ∗ λ ∗ − λ ^ 0 = γ ∗ − γ i ∗ λ ∗ ' γ ∗ − γ i ∗ λ ∗ p γ i ∗ γ i σ 2 = h v + k / 2 2 π N 2 exp ( − h 2 v ¯ S 2 + λ ∗ − λ ^ 0 ' γ i ' γ i ∗ λ ∗ − λ ^ 0 (27) v ¯ = N − K and N = v ¯ + K

Consider an informative prior created by set. v − 1 ¯ j = 1 c k j γ j ∗ γ j ∗

And letting c → 0 for j = 1 , 2

The posterior distribution of λ ∗ conditional on γ ∗ , h, Ω is given by: p λ ∗ γ ∗ h Ω α h N 2 exp ( − h 2 γ ∗ − λ ∗ γ ∗ ' Ω − 1 γ ∗ − λ ∗ γ ∗ + λ ∗ − λ 0 ∗ ' V − 1 ¯ γ ∗ − λ 0 ∗ × exp − h 2 Ω γ i ∗ − λ ∗ γ i ∗ ' γ i ∗ − λ ∗ γ i ∗ + λ ∗ − λ 0 ∗ ' λ ∗ − λ 0 ∗ v ¯ (28)

Solving the exponential part of the above equation, we will have: γ i ∗ − λ ∗ γ i ∗ ' γ i ∗ − λ ∗ γ i ∗ = γ ∗ 2 + λ ∗ γ ∗ 2 − 2 y i ∗ λ ∗ γ ∗ and γ ∗ − λ 0 ∗ ' λ ∗ − λ 0 ∗ = λ ∗ 2 + λ 0 ∗ 2 − 2 λ ∗ λ 0 ∗

Therefore, = exp − h 2 Ω ∑ i = 1 N y i ∗ 2 + λ ∗ γ ∗ 2 − 2 y i ∗ λ ∗ γ i ∗ + λ ∗ − λ 0 ∗ v ¯ 2

The additional term not involving λ ∗ is factored out to give: = exp − λ ∗ 2 2 v 2 ¯ + λ ∗ 2 λ 0 ∗ v 2 ¯ + λ ∗ ny ∗ Ω σ 2 − n λ ∗ 2 2 Ω σ 2 (29)

Factorization in terms of λ ∗ , the term in the exponential becomes: = − λ ∗ 2 σ 2 + 2 λλ ∗ ∗ 2 σ ∗ 2 σ ∗ 2 = 1 v ¯ + n Ω σ 2 − 1 and λ ∗ ∗ = σ ∗ 2 λ 0 2 v ¯ + ny ∗ Ω σ 2

So, the posterior density of λ ∗ conditioned on other parameter h, Ω , y ^∗ is a multivariate normal with mean λ ∗ and variance σ ∗ 2 .

That is, p λ ∗ h Ω y ∗ ∼ N λ ∗ ∗ σ ∗ 2

The posterior distribution of h conditional on λ ∗ , Ω , y ∗ is given by: P λ ∗ h Ω γ ∗ α h N 2 exp − h 2 γ ∗ − λ ∗ γ ∗ ' Ω − 1 γ ∗ − λ ∗ γ ∗ × h N + v − k 2 exp − hv 2 S 2 = h N + v − k 2 exp − h 2 Ω ∑ i N y i 2 ∗ + n λ ∗ γ ∗ 2 − 2 y ∗ n λ ∗ γ ∗ − hv 2 s 2 (30)

The posterior distribution of Ω ^*, conditional on y ^*, λ ^*, h, is given by: P Ω y ∗ λ ∗ h α P Ω × h N 2 exp − h 2 y ∗ − λ ∗ γ ∗ ′ Ω − 1 y ∗ − λ ∗ γ ∗ (31) P Ω y ∗ λ ∗ h α h N 2 exp − h 2 y ∗ − λ ∗ γ ∗ ′ Ω − 1 y ∗ − λ ∗ γ ∗ × | Ω ∗ | − ( β 0 + k + 1 ) / 2 exp tr R 0 − 1 β 0 − 1 / 2 (32)

The Gibbs sampling procedure used in this study involves generation of sequence of draws from the conditional posterior distribution of each parameter. ^{2 , 22 , 25}

(i)

Chose a starting or initial value, ϕ 0 for s = 1 , 2 , … , S

The right-hand side of ( 15) is proportional to the density function of an inverse Wishart distribution

Then, P Φ Y Ω ∼ IW q ΩΩ T R 0 − 1 n + ρ 0 (33)

•

At different functional forms of ³ heteroscedastic error structure with changes in sample size of 50, 100, 200 and 500. Hyper-parameter will be arbitrarily chosen for the simulation using Gibbs sampler an MCMC method. ^{6 , 22}

In order to evaluate the Bayesian model fit, we used the posterior predictive probability (PPP) procedure. ^{4 , 5 , 7 , 24} PPP = P ( f y λ ̂ i < f y rep λ ̂ i ≡ 1 m ∑ i = 1 m δ i (34)

After achieving convergence (after j iterations). λ ̂ ∗ j + 1 λ j + 1 Ω j + 1 can be regarded as observation from p(λ*, Ω| y) collect λ ∗ t Ω ∗ t t = j + 1 . … + T for statistical inference. λ ̂ = T − 1 ∑ t = 1 T λ t , Ω ̂ = T − 1 ∑ t = 1 T Ω t (35) gives Bayesian estimates of parameter and the latent variables. ^{10 , 17 , 23}

This gives the results for the latent and observed variables at various sample sizes for the four heteroscedastic error conditions considered.

Using the assumed values for the estimates which are λ 1 = 2.0, λ 2 = 3.0 and precision = 15.0.

The covariance matrix of ω was derived to be E ⋅ ηξT = ∏ 0 − 1 Γ with M at fixed values (0 or 1). The Bayesian estimates of SEM using the independent normal-gamma priors were derived for the two classes of SEM. Hyper-parameter was arbitrarily chosen for the simulation using Gibbs sampler a Markov chain Monte Carlo (MCMC) method since the joint posterior density does not have a tractable form. For the double logarithmic form, at 95% credible interval, when n=50, Posterior Mean, PM, and Precision, PR (2.011, 2.435, and 13.202), Posterior Standard Deviation PSD (0.035, 0.033, and 0.223) and when n=100, PM, and PR (2.022, 2.528, and 13.70), PSD (0.023, 0.025, and 0.251), when n=200, PM, and PR (2.052, 2.611, and 14.4), PSD (0.017, 0.018, and 0.255), when n=500, PM, and PR (2.010, 2.801, and 14.7), PSD (0.031, 0.021, and 0.258).

For the linear form, when n=50, PM, and PR (1.845, 2.779, and 13.95), PSD (0.240, 0.242, and 0.235). When n=100, PM, and PR (1.861, 2.811, and 14.22), PSD (0.328, 0.226, and 0.325), when n= 200, PM, and PR (1.956, 2.921, and 14.72), PSD (0.219, 0.217, and 0.212), and when n=500, PM, and PR (2.120, 3.122, and 14.95), PSD (0.211, 0.311, and 0.114).

For the linear-inverse form when n=50, PM, and PR (1.882, 2.742, and 14.95), PSD (0.040, 0.028, and 0.291). When n=100, PM, and PR (1.972, 2.835, and 14.65), PSD (0.024, 0.023, and 0.229). When n=200, PM, and PR (1.988, 2.901, and 14.45), PSD (0.017, 0.016, and 0.109), and when n=500, PM, and PR (2.021, 3.003, and 14.21), PSD (0.011, 0.015, and 0.105).

For the linear-absolute form, when n=50, PM, and PR (2.036, 2.824, and 14.500), PSD (0.032, 0.034, and 0.122), When n=100, PM, and PR (1.908, 2.903, and 13.92), PSD (0.022, 0.026, and 0.234). When n=200, PM, and PR (1.893, 2.809, and 13.85), PSD (0.017, 0.023, and 0.311), and when n=500, PM, and PR (1.806, 2.788, and 13.55), PSD (0.031, 0.035, and 0.433).

Examining different forms of heteroscedastic error structures in Bayesian structural equation modeling using informative priors, rather than assuming homogenous variance which is often a statistical fallacy in many studies. We compare the models’ posterior means and standard deviations in Tables 1, 2, 3 and 4. The differences are unlikely to impact substantive conclusions, but two of them are noteworthy.

First, the posterior means of the loadings ( λ 1 and λ 2 ) are somewhat smaller under different heteroscedastic condition with the informative priors as observed in Tables 6 and 7. Second, the factor variance γ ∗ is larger under our model with informative priors, likely because the informative prior placed more density on larger values of the posterior standard deviation. An evaluation of the model fit was based on the values of PPP as shown in Table 5 and it was observed that the linear form is the best with minimum PPP value as sample size increases. It was also revealed by the downward slope of the model as the sample size increases from 50 to 500 shown in Figure 1b when compared with Figure 1a, 2a and 2b.

Considering an improvement to maximum likelihood method, in Bayesian estimations, parameters are considered as random with informative prior distribution also known as the conjugate family of the posterior, once the data is simulated/collected, it is combined with prior distribution using Bayes theorem, next posterior distribution is calculated reflecting the prior knowledge and simulated data. ^{14 , 15 , 21} Joint posterior distribution is summarized using MCMC simulation techniques in terms of lower dimensional summary statistics as posterior mean and posterior standard deviations. ^{5 , 25} We observe that the structural and measurement equation obtained from this study are adequate and in general we could accept the proposed model.

All data underlying the results are available as part of the article and no additional source data are required.

Figshare: RCODE BSEM.docx. https://doi.org/10.6084/m9.figshare.19299851. ²⁶

Data are available under the terms of the Creative Commons Zero “No rights reserved” data waiver (CC0 1.0 Public domain dedication).