Fit Latent Data"

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Fit Latent Data

Enjoy this brief demonstration of the fit latent data module (i.e., a simple mediation model).

Also, please see Fit Observed Data.

First we simulate some data

# Number of observations
n <- 1000
# Coefficient for a path (x -> m)
a <- .75
# Coefficient for b path (m -> y)
b <- .80
# Coefficient for total effect (x -> y)
c <- .60
# Coefficient for indirect effect (x -> m -> y)
ab <- a * b
# Coefficient for direct effect (x -> y)
cd <- c - ab
# Compute x, m, y values
set.seed(100)
x <- rnorm(n)
m <- a * x + sqrt(1 - a^2) * rnorm(n)
eps <- 1 - (cd^2 + b^2 + 2*a * cd * b)
y <- cd * x + b * m + eps * rnorm(n)

data <- data.frame(y = y,
                   x = x,
                   m = m)

Next we run a standard mediation analysis using lavaan

model <- "
# direct effect
y ~ c*x
# mediator
m ~ a*x
y ~ b*m
# indirect effect (a*b)
ab := a*b
# total effect
cd := c + (a*b)"

fit <- lavaan::sem(model, data = data)
lavaan::summary(fit)

Now for the Bayesian model

bayesian.fit <- bfw::bfw(project.data = data,
                    latent = "x,m,y",
                    saved.steps = 50000,
                    latent.names = "Independent,Mediator,Dependent",
                    additional = "indirect <- xm * my , total <- xy + (xm * my)",
                    additional.names = "AB,C",
                    jags.model = "fit",
                    silent = TRUE)

round(bayesian.fit$summary.MCMC[,3:7],3)
#>                                       Mode   ESS  HDIlo HDIhi    n
#> beta[2,1]: Mediator vs. Independent  0.760 48988  0.720 0.799 1000
#> beta[3,1]: Dependent vs. Independent 0.024 13042 -0.012 0.058 1000
#> beta[3,2]: Dependent vs. Mediator    0.751 13230  0.715 0.786 1000
#> indirect[1]: AB                      0.571 21431  0.531 0.611 1000
#> total[1]: C                          0.591 49074  0.555 0.630 1000

Time for some noise

biased.sigma <-matrix(c(1,1,0,1,1,0,0,0,1),3,3)
set.seed(101)
noise <- MASS::mvrnorm(n=2,
                       mu=c(200, 300, 0),
                       Sigma=biased.sigma,
                       empirical=FALSE)
colnames(noise) <- c("y","x","m")
biased.data <- rbind(data,noise)

Rerun the lavaan model

biased.fit <- lavaan::sem(model, data = biased.data)
lavaan::summary(biased.fit)

Run the Bayesian model with robust estimates

biased.bfit <- bfw::bfw(project.data = data,
                    latent = "x,m,y",
                    saved.steps = 50000,
                    latent.names = "Independent,Mediator,Dependent",
                    additional = "indirect <- xm * my , total <- xy + (xm * my)",
                    additional.names = "AB,C",
                    jags.model = "fit",
                    run.robust = TRUE,
                    jags.seed = 101,
                    silent = TRUE)

round(biased.bfit$summary.MCMC[,3:7],3)
#>                                       Mode   ESS  HDIlo HDIhi    n
#> beta[2,1]: Mediator vs. Independent  0.763 31178  0.721 0.799 1000
#> beta[3,1]: Dependent vs. Independent 0.022  7724 -0.014 0.057 1000
#> beta[3,2]: Dependent vs. Mediator    0.751  7772  0.714 0.786 1000
#> indirect[1]: AB                      0.572 12913  0.531 0.610 1000
#> total[1]: C                          0.590 31362  0.557 0.631 1000


Try the bfw package in your browser

Any scripts or data that you put into this service are public.

bfw documentation built on Jan. 11, 2020, 9:29 a.m.