In oeysan/bfw: Bayesian Framework for Computational Modeling
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
Fit Observed Data
Enjoy this brief demonstration of the fit observed data module (i.e., a simple mediation model).
Also, please see Fit Latent 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)
Then we create some observed data based on the latent variables
set.seed(102)
## create 5 y variables, 5 x variables and 5 m variables with jittered data
jitter.data <- lapply(1:5, function (i) {
apply(data , 2 , function (x) jitter(x , amount=1))
})
# Create data frame with jittered observed variables
jitter.data <- data.frame( sapply(jitter.data, "[" ,,1) , # First column (y)
sapply(jitter.data, "[" ,,2) , # Second column (x)
sapply(jitter.data, "[" ,,3) # Third column (m)
)
# Add column names
colnames(jitter.data) <- c( paste0("y",1:5) , paste0("x",1:5) , paste0("m",1:5))
Next we run a standard observed mediation analysis using lavaan
model <- '
y =~ y1+y2+y3+y4+y5
x =~ x1+x2+x3+x4+x5
m =~ m1+m2+m3+m4+m5
# 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)
'
jitter.fit <- lavaan::sem(model, data = jitter.data)
lavaan::summary(jitter.fit)
Run the Bayesian model on observed data
bayesian.jitter.fit <- bfw(project.data = jitter.data,
observed = "x1,x2,x3,x4,x5,
m1,m2,m3,m4,m5,
y1,y2,y3,y4,y5",
latent.names = "Independent,Mediator,Dependent",
additional = "indirect <- xm * my , total <- xy + (xm * my)",
additional.names = "AB,C",
jags.model = "fit",
jags.seed = 102,
silent = TRUE)
round(bayesian.jitter.fit$summary.MCMC[,3:7],3)
#> Mode ESS HDIlo HDIhi n
#> lam[1]: X1 0.999 0 1.000 1.000 1000
#> lam[2]: X2 1.019 1309 0.966 1.067 1000
#> lam[3]: X3 0.988 1420 0.941 1.042 1000
#> lam[4]: X4 0.998 1440 0.946 1.048 1000
#> lam[5]: X5 0.990 1402 0.932 1.035 1000
#> lam[6]: M1 0.999 0 1.000 1.000 1000
#> lam[7]: M2 0.964 1164 0.912 1.013 1000
#> lam[8]: M3 0.970 1049 0.927 1.026 1000
#> lam[9]: M4 1.013 1075 0.959 1.062 1000
#> lam[10]: M5 0.969 1117 0.921 1.022 1000
#> lam[11]: Y1 0.999 0 1.000 1.000 1000
#> lam[12]: Y2 1.004 1136 0.939 1.059 1000
#> lam[13]: Y3 0.994 1139 0.933 1.050 1000
#> lam[14]: Y4 0.941 1270 0.883 0.999 1000
#> lam[15]: Y5 1.023 1149 0.954 1.074 1000
#> error[1]: X1 0.304 4691 0.273 0.340 1000
#> error[2]: X2 0.321 4464 0.285 0.355 1000
#> error[3]: X3 0.327 4637 0.296 0.366 1000
#> error[4]: X4 0.343 5200 0.305 0.379 1000
#> error[5]: X5 0.346 4795 0.309 0.383 1000
#> error[6]: M1 0.352 5155 0.315 0.388 1000
#> error[7]: M2 0.329 5256 0.296 0.365 1000
#> error[8]: M3 0.319 4979 0.288 0.355 1000
#> error[9]: M4 0.321 4903 0.291 0.358 1000
#> error[10]: M5 0.322 5080 0.288 0.356 1000
#> error[11]: Y1 0.300 4509 0.272 0.337 1000
#> error[12]: Y2 0.333 4645 0.298 0.367 1000
#> error[13]: Y3 0.317 5032 0.288 0.356 1000
#> error[14]: Y4 0.350 5481 0.317 0.388 1000
#> error[15]: Y5 0.329 4929 0.294 0.364 1000
#> cov[1,1]: Independent vs. Independent 1.044 884 0.968 1.114 1000
#> cov[2,1]: Mediator vs. Independent 0.822 949 0.781 0.872 1000
#> cov[3,1]: Dependent vs. Independent 0.609 936 0.574 0.647 1000
#> cov[1,2]: Independent vs. Mediator 0.822 949 0.781 0.872 1000
#> cov[2,2]: Mediator vs. Mediator 1.071 669 0.994 1.151 1000
#> cov[3,2]: Dependent vs. Mediator 0.798 776 0.756 0.848 1000
#> cov[1,3]: Independent vs. Dependent 0.609 936 0.574 0.647 1000
#> cov[2,3]: Mediator vs. Dependent 0.798 776 0.756 0.848 1000
#> cov[3,3]: Dependent vs. Dependent 0.759 683 0.701 0.827 1000
#> beta[2,1]: Mediator vs. Independent 0.788 1308 0.735 0.850 1000
#> beta[3,1]: Dependent vs. Independent -0.014 1008 -0.072 0.043 1000
#> beta[3,2]: Dependent vs. Mediator 0.760 630 0.696 0.828 1000
#> indirect[1]: AB 0.593 917 0.535 0.662 1000
#> total[1]: C 0.587 1713 0.534 0.640 1000
#> Fit (Predicted) 80.578 5928 69.388 97.268 1000
#> Fit (Simulated) 119.033 10000 92.518 152.408 1000
#> Fit (Discrepancy) -38.274 9228 -71.546 -5.316 1000
#> PPP 0.991 20 0.988 0.992 1000
Let's add some noise to the observed data
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)
jitter.noise <- noise[,rep(1:3,each=5)]
colnames(jitter.noise) <- c( paste0("y",1:5) , paste0("x",1:5) , paste0("m",1:5))
biased.jitter.data <- rbind(jitter.data, jitter.noise)
Run the lavaan model on biased data
biased.jitter.fit <- lavaan::sem(model, data = biased.jitter.data)
lavaan::summary(biased.jitter.fit)
Finally, run the Bayesian model with robust estimates on biased data
bayesian.biased.jitter.fit <- bfw(project.data = biased.jitter.data,
observed = "x1,x2,x3,x4,x5,
m1,m2,m3,m4,m5,
y1,y2,y3,y4,y5",
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 = 103,
silent = TRUE)
round(bayesian.biased.jitter.fit$summary.MCMC[,3:7],3)
#> Mode ESS HDIlo HDIhi n
#> lam[1]: X1 0.999 0 1.000 1.000 1002
#> lam[2]: X2 1.011 813 0.961 1.068 1002
#> lam[3]: X3 0.983 846 0.933 1.037 1002
#> lam[4]: X4 0.987 833 0.938 1.048 1002
#> lam[5]: X5 0.973 906 0.925 1.032 1002
#> lam[6]: M1 0.999 0 1.000 1.000 1002
#> lam[7]: M2 0.961 633 0.911 1.014 1002
#> lam[8]: M3 0.977 624 0.925 1.026 1002
#> lam[9]: M4 1.014 588 0.958 1.060 1002
#> lam[10]: M5 0.969 630 0.921 1.022 1002
#> lam[11]: Y1 0.999 0 1.000 1.000 1002
#> lam[12]: Y2 0.998 765 0.940 1.062 1002
#> lam[13]: Y3 1.012 816 0.943 1.068 1002
#> lam[14]: Y4 0.938 936 0.876 0.999 1002
#> lam[15]: Y5 1.023 717 0.962 1.090 1002
#> error[1]: X1 0.223 1986 0.186 0.255 1002
#> error[2]: X2 0.230 2228 0.200 0.275 1002
#> error[3]: X3 0.249 2177 0.214 0.287 1002
#> error[4]: X4 0.263 2074 0.223 0.299 1002
#> error[5]: X5 0.271 2583 0.232 0.308 1002
#> error[6]: M1 0.343 3209 0.309 0.381 1002
#> error[7]: M2 0.324 2957 0.289 0.359 1002
#> error[8]: M3 0.312 2928 0.282 0.352 1002
#> error[9]: M4 0.315 2939 0.283 0.352 1002
#> error[10]: M5 0.319 3258 0.284 0.352 1002
#> error[11]: Y1 0.221 2172 0.192 0.260 1002
#> error[12]: Y2 0.259 2553 0.226 0.297 1002
#> error[13]: Y3 0.248 2543 0.209 0.279 1002
#> error[14]: Y4 0.284 2934 0.247 0.322 1002
#> error[15]: Y5 0.250 2396 0.215 0.288 1002
#> cov[1,1]: Independent vs. Independent 1.051 543 0.974 1.126 1002
#> cov[2,1]: Mediator vs. Independent 0.825 597 0.782 0.872 1002
#> cov[3,1]: Dependent vs. Independent 0.604 603 0.566 0.643 1002
#> cov[1,2]: Independent vs. Mediator 0.825 597 0.782 0.872 1002
#> cov[2,2]: Mediator vs. Mediator 1.070 377 0.997 1.159 1002
#> cov[3,2]: Dependent vs. Mediator 0.797 541 0.746 0.838 1002
#> cov[1,3]: Independent vs. Dependent 0.604 603 0.566 0.643 1002
#> cov[2,3]: Mediator vs. Dependent 0.797 541 0.746 0.838 1002
#> cov[3,3]: Dependent vs. Dependent 0.761 458 0.701 0.832 1002
#> beta[2,1]: Mediator vs. Independent 0.787 715 0.732 0.847 1002
#> beta[3,1]: Dependent vs. Independent -0.024 621 -0.080 0.041 1002
#> beta[3,2]: Dependent vs. Mediator 0.750 420 0.684 0.825 1002
#> indirect[1]: AB 0.588 604 0.528 0.660 1002
#> total[1]: C 0.572 1358 0.523 0.629 1002
oeysan/bfw documentation built on March 27, 2022, 7:41 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
Fit Observed Data
Enjoy this brief demonstration of the fit observed data module (i.e., a simple mediation model).
Also, please see Fit Latent 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)
Then we create some observed data based on the latent variables
set.seed(102) ## create 5 y variables, 5 x variables and 5 m variables with jittered data jitter.data <- lapply(1:5, function (i) { apply(data , 2 , function (x) jitter(x , amount=1)) }) # Create data frame with jittered observed variables jitter.data <- data.frame( sapply(jitter.data, "[" ,,1) , # First column (y) sapply(jitter.data, "[" ,,2) , # Second column (x) sapply(jitter.data, "[" ,,3) # Third column (m) ) # Add column names colnames(jitter.data) <- c( paste0("y",1:5) , paste0("x",1:5) , paste0("m",1:5))
Next we run a standard observed mediation analysis using lavaan
model <- ' y =~ y1+y2+y3+y4+y5 x =~ x1+x2+x3+x4+x5 m =~ m1+m2+m3+m4+m5 # 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) ' jitter.fit <- lavaan::sem(model, data = jitter.data) lavaan::summary(jitter.fit)
Run the Bayesian model on observed data
bayesian.jitter.fit <- bfw(project.data = jitter.data, observed = "x1,x2,x3,x4,x5, m1,m2,m3,m4,m5, y1,y2,y3,y4,y5", latent.names = "Independent,Mediator,Dependent", additional = "indirect <- xm * my , total <- xy + (xm * my)", additional.names = "AB,C", jags.model = "fit", jags.seed = 102, silent = TRUE) round(bayesian.jitter.fit$summary.MCMC[,3:7],3) #> Mode ESS HDIlo HDIhi n #> lam[1]: X1 0.999 0 1.000 1.000 1000 #> lam[2]: X2 1.019 1309 0.966 1.067 1000 #> lam[3]: X3 0.988 1420 0.941 1.042 1000 #> lam[4]: X4 0.998 1440 0.946 1.048 1000 #> lam[5]: X5 0.990 1402 0.932 1.035 1000 #> lam[6]: M1 0.999 0 1.000 1.000 1000 #> lam[7]: M2 0.964 1164 0.912 1.013 1000 #> lam[8]: M3 0.970 1049 0.927 1.026 1000 #> lam[9]: M4 1.013 1075 0.959 1.062 1000 #> lam[10]: M5 0.969 1117 0.921 1.022 1000 #> lam[11]: Y1 0.999 0 1.000 1.000 1000 #> lam[12]: Y2 1.004 1136 0.939 1.059 1000 #> lam[13]: Y3 0.994 1139 0.933 1.050 1000 #> lam[14]: Y4 0.941 1270 0.883 0.999 1000 #> lam[15]: Y5 1.023 1149 0.954 1.074 1000 #> error[1]: X1 0.304 4691 0.273 0.340 1000 #> error[2]: X2 0.321 4464 0.285 0.355 1000 #> error[3]: X3 0.327 4637 0.296 0.366 1000 #> error[4]: X4 0.343 5200 0.305 0.379 1000 #> error[5]: X5 0.346 4795 0.309 0.383 1000 #> error[6]: M1 0.352 5155 0.315 0.388 1000 #> error[7]: M2 0.329 5256 0.296 0.365 1000 #> error[8]: M3 0.319 4979 0.288 0.355 1000 #> error[9]: M4 0.321 4903 0.291 0.358 1000 #> error[10]: M5 0.322 5080 0.288 0.356 1000 #> error[11]: Y1 0.300 4509 0.272 0.337 1000 #> error[12]: Y2 0.333 4645 0.298 0.367 1000 #> error[13]: Y3 0.317 5032 0.288 0.356 1000 #> error[14]: Y4 0.350 5481 0.317 0.388 1000 #> error[15]: Y5 0.329 4929 0.294 0.364 1000 #> cov[1,1]: Independent vs. Independent 1.044 884 0.968 1.114 1000 #> cov[2,1]: Mediator vs. Independent 0.822 949 0.781 0.872 1000 #> cov[3,1]: Dependent vs. Independent 0.609 936 0.574 0.647 1000 #> cov[1,2]: Independent vs. Mediator 0.822 949 0.781 0.872 1000 #> cov[2,2]: Mediator vs. Mediator 1.071 669 0.994 1.151 1000 #> cov[3,2]: Dependent vs. Mediator 0.798 776 0.756 0.848 1000 #> cov[1,3]: Independent vs. Dependent 0.609 936 0.574 0.647 1000 #> cov[2,3]: Mediator vs. Dependent 0.798 776 0.756 0.848 1000 #> cov[3,3]: Dependent vs. Dependent 0.759 683 0.701 0.827 1000 #> beta[2,1]: Mediator vs. Independent 0.788 1308 0.735 0.850 1000 #> beta[3,1]: Dependent vs. Independent -0.014 1008 -0.072 0.043 1000 #> beta[3,2]: Dependent vs. Mediator 0.760 630 0.696 0.828 1000 #> indirect[1]: AB 0.593 917 0.535 0.662 1000 #> total[1]: C 0.587 1713 0.534 0.640 1000 #> Fit (Predicted) 80.578 5928 69.388 97.268 1000 #> Fit (Simulated) 119.033 10000 92.518 152.408 1000 #> Fit (Discrepancy) -38.274 9228 -71.546 -5.316 1000 #> PPP 0.991 20 0.988 0.992 1000
Let's add some noise to the observed data
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) jitter.noise <- noise[,rep(1:3,each=5)] colnames(jitter.noise) <- c( paste0("y",1:5) , paste0("x",1:5) , paste0("m",1:5)) biased.jitter.data <- rbind(jitter.data, jitter.noise)
Run the lavaan model on biased data
biased.jitter.fit <- lavaan::sem(model, data = biased.jitter.data) lavaan::summary(biased.jitter.fit)
Finally, run the Bayesian model with robust estimates on biased data
bayesian.biased.jitter.fit <- bfw(project.data = biased.jitter.data, observed = "x1,x2,x3,x4,x5, m1,m2,m3,m4,m5, y1,y2,y3,y4,y5", 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 = 103, silent = TRUE) round(bayesian.biased.jitter.fit$summary.MCMC[,3:7],3) #> Mode ESS HDIlo HDIhi n #> lam[1]: X1 0.999 0 1.000 1.000 1002 #> lam[2]: X2 1.011 813 0.961 1.068 1002 #> lam[3]: X3 0.983 846 0.933 1.037 1002 #> lam[4]: X4 0.987 833 0.938 1.048 1002 #> lam[5]: X5 0.973 906 0.925 1.032 1002 #> lam[6]: M1 0.999 0 1.000 1.000 1002 #> lam[7]: M2 0.961 633 0.911 1.014 1002 #> lam[8]: M3 0.977 624 0.925 1.026 1002 #> lam[9]: M4 1.014 588 0.958 1.060 1002 #> lam[10]: M5 0.969 630 0.921 1.022 1002 #> lam[11]: Y1 0.999 0 1.000 1.000 1002 #> lam[12]: Y2 0.998 765 0.940 1.062 1002 #> lam[13]: Y3 1.012 816 0.943 1.068 1002 #> lam[14]: Y4 0.938 936 0.876 0.999 1002 #> lam[15]: Y5 1.023 717 0.962 1.090 1002 #> error[1]: X1 0.223 1986 0.186 0.255 1002 #> error[2]: X2 0.230 2228 0.200 0.275 1002 #> error[3]: X3 0.249 2177 0.214 0.287 1002 #> error[4]: X4 0.263 2074 0.223 0.299 1002 #> error[5]: X5 0.271 2583 0.232 0.308 1002 #> error[6]: M1 0.343 3209 0.309 0.381 1002 #> error[7]: M2 0.324 2957 0.289 0.359 1002 #> error[8]: M3 0.312 2928 0.282 0.352 1002 #> error[9]: M4 0.315 2939 0.283 0.352 1002 #> error[10]: M5 0.319 3258 0.284 0.352 1002 #> error[11]: Y1 0.221 2172 0.192 0.260 1002 #> error[12]: Y2 0.259 2553 0.226 0.297 1002 #> error[13]: Y3 0.248 2543 0.209 0.279 1002 #> error[14]: Y4 0.284 2934 0.247 0.322 1002 #> error[15]: Y5 0.250 2396 0.215 0.288 1002 #> cov[1,1]: Independent vs. Independent 1.051 543 0.974 1.126 1002 #> cov[2,1]: Mediator vs. Independent 0.825 597 0.782 0.872 1002 #> cov[3,1]: Dependent vs. Independent 0.604 603 0.566 0.643 1002 #> cov[1,2]: Independent vs. Mediator 0.825 597 0.782 0.872 1002 #> cov[2,2]: Mediator vs. Mediator 1.070 377 0.997 1.159 1002 #> cov[3,2]: Dependent vs. Mediator 0.797 541 0.746 0.838 1002 #> cov[1,3]: Independent vs. Dependent 0.604 603 0.566 0.643 1002 #> cov[2,3]: Mediator vs. Dependent 0.797 541 0.746 0.838 1002 #> cov[3,3]: Dependent vs. Dependent 0.761 458 0.701 0.832 1002 #> beta[2,1]: Mediator vs. Independent 0.787 715 0.732 0.847 1002 #> beta[3,1]: Dependent vs. Independent -0.024 621 -0.080 0.041 1002 #> beta[3,2]: Dependent vs. Mediator 0.750 420 0.684 0.825 1002 #> indirect[1]: AB 0.588 604 0.528 0.660 1002 #> total[1]: C 0.572 1358 0.523 0.629 1002
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.