Bayesian methods for estimating structural equation models

1. Basic concepts

Let \(M\) be an arbitrary SEM with a vector of unknown parameters \(\mathbf{\theta}\),
Let \(\mathbf{Y} = (\mathbf{y}_1, \dots , \mathbf{y}_n)\) be the observed data set of raw observations with a sample size \(n\).
In a Bayesian approach, \(\mathbf{\theta}\) is considered to be random with a distribution (called a prior distribution) and an associated (prior) density function, say, \(p(\mathbf{\theta}|M)\). For simplicity, we use \(p(\mathbf{\theta})\) to denote \(p(\mathbf{\theta}|M)\).
Bayesian estimation is based on the observed data \(\mathbf{Y}\) and the prior distribution of \(\mathbf{\theta}\).
Let \(p(\mathbf{Y},\mathbf{\theta}|M)\) be the probability density function of the joint distribution of \(\mathbf{Y}\) and \(\mathbf{\theta}\) under \(M\).
The behavior of \(\mathbf{\theta}\) under the given data \(\mathbf{Y}\) is fully described by the posterior density function of \(\mathbf{\theta}\) defined by. Let \(p(\mathbf{\theta},|\mathbf{Y}, M)\).
Based on a wellknown identity in probability, we have \[ log p(\mathbf{\theta},|\mathbf{Y}, M)=log p(\mathbf{Y}|\mathbf{\theta}, M) +log p(\mathbf{\theta})+ constant. \]

2. Specification of priors

Consider

Measurement equation: \[\mathbf{y}_i=\mathbf{\mu}+\mathbf{\Lambda}\mathbf{\omega}_i +\mathbf{\epsilon}_i\],

Structural equation: \[\mathbf{\eta}_i=\mathbf{B}\mathbf{d}_i+\mathbf{\Pi}\mathbf{\eta}_i+\mathbf{\Gamma}\mathbf{F}(\mathbf{\xi}_i)+ \mathbf{\delta}_i\] Consider the following priors specification from authors

Fig 1. Priors specification

where hyperparameters are specified following

Fig 2. Hyperparameters

Of course, also non-conjugate priors can be proposed for the Bayesian estimation of SEM. But, since that posterior priors

Theoretically, it could be obtained via integration. For most situations, the integration does not have a closed form. However, if we can simulate a sufficiently large number of observations from posterior distributions, we can approximate the mean and other useful statistics through the simulated observations. Hence, to solve the problem, it suffices to develop efficient and dependable methods for drawing observations from the posterior distribution. For most nonstandard SEMs, the posterior distribution is complicated.

More specifically, instead of working on the intractable posterior density \(p(\mathbf{\theta},|\mathbf{Y})\), we will work on \(p(\mathbf{\theta},\mathbf{\Omega}|\mathbf{Y})\), where \(\mathbf{\Omega}\) is the set of latent variables in the model in a augmented version of the model. For most cases, \(p(\mathbf{\theta},\mathbf{\Omega}|\mathbf{Y})\) is still not in closed form and it is difficult to deal with it directly. However, on the basis of the complete data set \((\mathbf{\Omega}|\mathbf{Y})\), the conditional distribution \(p(\mathbf{\theta}|\mathbf{\Omega}, \mathbf{Y})\) is usually standard, and the conditional distribution p(|,) can also be derived from the definition of the model without much difficulty. As a result, we can apply some MCMC methods.

3. Bayesian estimation via WinBUGS

if(!require(mvtnorm)) install.packages("mvtnorm");
if(!require(R2WinBUGS)) install.packages("R2WinBUGS");
library(mvtnorm)   #Load mvtnorm package
library(R2WinBUGS) #Load R2WinBUGS package

N=500                         #Sample size
BZ=numeric(N)                 #Fixed covariate in structural equation
XI=matrix(NA, nrow=N, ncol=2) #Explanatory latent variables
Eta=numeric(N)                #Outcome latent variables
Y=matrix(NA, nrow=N, ncol=8)  #Observed variables

#The covariance matrix of xi
phi=matrix(c(1, 0.3, 0.3, 1), nrow=2)

#Estimates and standard error estimates
Eu=matrix(NA, nrow=100, ncol=10);    SEu=matrix(NA, nrow=100, ncol=10)
Elam=matrix(NA, nrow=100, ncol=7);   SElam=matrix(NA, nrow=100, ncol=7)
Eb=numeric(100);                     SEb=numeric(100)
Egam=matrix(NA, nrow=100, ncol=5);   SEgam=matrix(NA, nrow=100, ncol=5)
Esgm=matrix(NA, nrow=100, ncol=10);  SEsgm=matrix(NA, nrow=100, ncol=10)
Esgd=numeric(100);                   SEsgd=numeric(100)
Ephx=matrix(NA, nrow=100, ncol=3);   SEphx=matrix(NA, nrow=100, ncol=3)

R=matrix(c(1.0, 0.3, 0.3, 1.0), nrow=2)

parameters=c("u", "lam", "b", "gam", "sgm", "sgd", "phx")

init1=list(u=rep(0,10), lam=rep(0,7), b=0, gam=rep(0,5), psi=rep(1,10),
           psd=1, phi=matrix(c(1, 0, 0, 1), nrow=2))

init2=list(u=rep(1,10), lam=rep(1,7), b=1, gam=rep(1,5), psi=rep(2,10),
           psd=2, phi=matrix(c(2, 0, 0, 2), nrow=2))

inits=list(init1, init2)

eps=numeric(10)

for (t in 1:100) {
    #Generate Data
    for (i in 1:N) {
        BZ[i]=rt(1, 5)

        XI[i,]=rmvnorm(1, c(0,0), phi)

        delta=rnorm(1, 0, sqrt(0.36))
        Eta[i]=0.5*BZ[i]+0.4*XI[i,1]+0.4*XI[i,2]+0.3*XI[i,1]*XI[i,2]
               +0.2*XI[i,1]*XI[i,1]+0.5*XI[i,2]*XI[i,2]+delta

        eps[1:3]=rnorm(3, 0, sqrt(0.3))
        eps[4:7]=rnorm(4, 0, sqrt(0.5))
        eps[8:10]=rnorm(3, 0, sqrt(0.4))
        Y[i,1]=Eta[i]+eps[1]
        Y[i,2]=0.9*Eta[i]+eps[2]
        Y[i,3]=0.7*Eta[i]+eps[3]
        Y[i,4]=XI[i,1]+eps[4]
        Y[i,5]=0.9*XI[i,1]+eps[5]
        Y[i,6]=0.7*XI[i,1]+eps[6]
        Y[i,7]=0.5*XI[i,1]+eps[7]
        Y[i,8]=XI[i,2]+eps[8]
        Y[i,9]=0.9*XI[i,2]+eps[9]
        Y[i,10]=0.7*XI[i,2]+eps[10]
    }

    #Run WinBUGS
    data=list(N=500, zero=c(0,0), z=BZ, R=R, y=Y)

    model<-bugs(data,inits,parameters,
                model.file="C:/Simulation/model.txt",
                n.chains=2,n.iter=10000,n.burnin=4000,n.thin=1,
                bugs.directory="c:/Program Files/WinBUGS14/",
                working.directory="C:/Simulation/")

    #Save Estimates
    Eu[t,]=model$mean$u;               SEu[t,]=model$sd$u
    Elam[t,]=model$mean$lam;           SElam[t,]=model$sd$lam
    Eb[t]=model$mean$b;                SEb[t]=model$sd$b
    Egam[t,]=model$mean$gam;           SEgam[t,]=model$sd$gam
    Esgm[t,]=model$mean$sgm;           SEsgm[t,]=model$sd$sgm
    Esgd[t]=model$mean$sgd;            SEsgd[t]=model$sd$sgd
    Ephx[t,1]=model$mean$phx[1,1];     SEphx[t,1]=model$sd$phx[1,1]
    Ephx[t,2]=model$mean$phx[1,2];     SEphx[t,2]=model$sd$phx[1,2]
    Ephx[t,3]=model$mean$phx[2,2];     SEphx[t,3]=model$sd$phx[2,2]
}

4. Results

alt text here

7. Observations

Data and software can be obtained in