Skip to contents

Inference-Example

Inference-Example.Rmd

The following shows a basic comparison of the inference feature provided by the Bayesian Network translation of bnSEM and a rejection sampling approach.

The objective is to investigate the distribution of a latent variable η\eta given different values of one of the indicators (x2x_2) of a second latent variable ξ\xi.

We want to fit the following model

model <- ' 
   xi =~ x1 + x2 + x3
   eta =~ y1 + y2 + y3
   eta ~ xi
   
   x2~~x3
'

To this end, we simulate some data:

# Simulation function:
simulate_data <- function(n){
  xi <- rnorm(n)
  eta <- .4*xi + rnorm(n)
  
  unobs <- rnorm(n)
  
  x1 <- 1*xi + rnorm(n,0,.2)
  x2 <- .7*xi + .3*unobs + rnorm(n,0,.2)
  x3 <- .5*xi + .2*unobs + rnorm(n,0,.2)
  
  y1 <- 1*eta + rnorm(n,0,.2)
  y2 <- .7*eta + rnorm(n,0,.2)
  y3 <- .5*eta + rnorm(n,0,.2)
  
  return(data.frame(xi, eta,
                    x1, x2, x3, y1, y2, y3))
}

data_set <- simulate_data(n = 100)

First, let’s use bnSEM to find the distribution of η\eta for x2∈(−2,−1)x_2 \in (-2, -1) and x2∈(1,3)x_2 \in (1,3):

# Now, we check if we can get the same distributions using a Bayesian Network:
fit <- mxsem(model, 
             data = data_set[,c("x1", "x2", "x3",
                                "y1", "y2", "y3")]) |>
  OpenMx::mxTryHard()
#> Running untitled1 with 20 parameters
#> 
#> Beginning initial fit attempt
#> Running untitled1 with 20 parameters
#> 
#>  Lowest minimum so far:  587.039224234481
#> 
#> Solution found
#> 
#>  Solution found!  Final fit=587.03922 (started at 2370.7554)  (1 attempt(s): 1 valid, 0 errors)
#>  Start values from best fit:
#> 0.58957322846888,0.437456148851343,0.362483373166715,0.70875447376711,0.484343167810718,0.13069588666507,0.0021926799799419,0.0749159820500079,0.0839271246208961,0.0575835519258439,0.032384521260863,0.046399855919389,0.84888770062111,0.921460926590316,0.120593804111261,0.0831610913840666,0.0608359960286526,-0.101313354961244,-0.0287944446546238,-0.0169745721886592

network <- bnSEM::bnSEM(mx_model = fit)
#> Found covariances in your model. The model will be translated to an equivalent model with phantom variables. The model with phantom variables will be returned in the internal list as 'internal_model'
#> Refitting the model with phantom variables. The fit will be the same, but the variances of the residuals will change.
#> Running untitled1 with 20 parameters
#> 
#> Beginning initial fit attempt
#> Running untitled1 with 20 parameters
#> 
#>  Lowest minimum so far:  589.436688213814
#>  Not all eigenvalues of the Hessian are positive: 68613.729749084, 40466.8107927976, 32314.9390672524, 16926.9063518055, 6904.53259090789, 6837.24043558328, 5143.62366210038, 4467.90310739625, 4125.23404759069, 3189.53521103073, 2259.4995508706, 1749.57121932741, 1645.13701356174, 871.302694638804, 180.661107400727, 141.535541064496, 109.888505320817, 98.0519312527793, 94.5109757487781, -101.117062150407
#> 
#> Beginning fit attempt 1 of at maximum 10 extra tries
#> Running untitled1 with 20 parameters
#> 
#>  Lowest minimum so far:  589.436688213812
#>  Not all eigenvalues of the Hessian are positive: 68613.7619321129, 40466.8063101584, 32314.9445017495, 16926.9210520897, 6904.53600570646, 6837.25929327516, 5143.62787655884, 4467.91730240462, 4125.23474395258, 3189.53764062725, 2259.49671686691, 1749.57484232925, 1645.14450331589, 871.304368966483, 180.662903773284, 141.543128954122, 109.887871547636, 98.055019407652, 94.5095611155132, -101.113711701767
#> 
#> Beginning fit attempt 2 of at maximum 10 extra tries
#> Running untitled1 with 20 parameters
#>  Not all eigenvalues of the Hessian are positive: 68613.7619321129, 40466.8063101584, 32314.9445017495, 16926.9210520897, 6904.53600570646, 6837.25929327516, 5143.62787655884, 4467.91730240462, 4125.23474395258, 3189.53764062725, 2259.49671686691, 1749.57484232925, 1645.14450331589, 871.304368966483, 180.662903773284, 141.543128954122, 109.887871547636, 98.055019407652, 94.5095611155132, -101.113711701767
#> 
#> Beginning fit attempt 3 of at maximum 10 extra tries
#> Running untitled1 with 20 parameters
#> 
#>  Lowest minimum so far:  587.039224234405
#> 
#> Solution found
#> 
#>  Solution found!  Final fit=587.03922 (started at 658.39424)  (4 attempt(s): 4 valid, 0 errors)
#>  Start values from best fit:
#> 0.589572540488107,0.437455695141567,0.362483898076655,0.708754632039264,0.484343246915199,-0.273708341974104,0.0557799788927127,0.00219166365995217,0.00901109139671079,0.0575835454247901,0.0323845658566736,0.0463998241209689,0.848888155806502,0.921461151361845,0.120593906200774,0.0831613432801809,0.0608360574459778,-0.101313138037528,-0.0287943593522188,-0.0169742017776087

dist_1 <- bnlearn::cpdist(fitted = network$bayes_net,
                          node = "eta",
                          evidence = (x2 > -2 & x2 < -1))
dist_2 <- bnlearn::cpdist(fitted = network$bayes_net,
                          node = "eta",
                          evidence = (x2 > 1 & x2 < 3))

We want to compare the distributions to rejection sampling, where we simply simulate large samples and only keep those subjects that meet our criterion:

target_n <- 10000
rejection_sampling <- c()
n <- 0
max_it <- 10
it <- 1
while((n < target_n) & (it < max_it)){
  data_ <- simulate_data(n = 1000000)
  rejection_sampling_1 <- rbind(rejection_sampling,
                                data_[data_$x2 > -2 & data_$x2 < -1,])
  rejection_sampling_2 <- rbind(rejection_sampling,
                                data_[data_$x2 > 1 & data_$x2 < 3,])
  n <- min(nrow(rejection_sampling_1), 
           nrow(rejection_sampling_2))
  it <- it+1
}

Comparing bnSEM and rejection sampling, we find the following distributions for η\eta:

combined <- rbind(data.frame(Type = "BN inference", 
                             Set = paste0("x2 in (", -2, ", ", -1, ")"), 
                             dist_1),
                  data.frame(Type = "BN inference", 
                             Set = paste0("x2 in (", 1, ", ", 3, ")"), 
                             dist_2),
                  data.frame(Type = "rejection sampling",
                             Set = paste0("x2 in (", -2, ", ", -1, ")"), 
                             eta = rejection_sampling_1[, "eta"]),
                  data.frame(Type = "rejection sampling",
                             Set = paste0("x2 in (", 1, ", ", 3, ")"), 
                             eta = rejection_sampling_2[, "eta"])
)

combined$Type <- as.factor(combined$Type)
combined$Set <- as.factor(combined$Set)

ggplot(combined, 
       aes(x = eta, color = Type, lty = Set)) +
  geom_density() +
  xlab(expression(eta)) +
  theme_classic()