Comparison AIC in linear regressions and path analysis

linear_and_path.R

# https://github.com/OhkuboYusaku/blog/tree/main/causal

library(lavaan)
library(parallel)

N.sim <- 10000
N.sample <- 100

sig.y <- 10

model1 <- '
  y ~ x
'

model2 <- '
  y ~ x + z1
'

model3 <- '
  y ~ x + z2
'

model4 <- '
  y  ~ x + z1
  x  ~ z2
  z1 ~ z2
'

# Data

set.seed(123)
data <- lapply(1:N.sim, function(i) {
  z2 <- rnorm(N.sample, 0, 10) # covariate
  z1 <- 1.2 * z2 + rnorm(N.sample, 0, 1) # covariate
  w <- -0.5 * z1 + rnorm(N.sample, 0, 1)
  x <-  0.5 * z2 + rnorm(N.sample, 0, 1) # target variable
  y <-  0.1 * w - 0.5 * x + 0.5 * z1 + rnorm(N.sample, 0, sig.y) # outcome

  data.frame(x = x, y = y, z1 = z1, z2 = z2)})

# Simulation

result <- mcmapply(function(d) {
  fit1 <- sem(model1, d)
  fit2 <- sem(model2, d)
  fit3 <- sem(model3, d)
  fit4 <- sem(model4, d)
  coef4 <- coef(fit4)
  c(AIC(fit1), AIC(fit2), AIC(fit3), AIC(fit4),
    coef4["y~x"], coef4["y~z1"])
},
data, mc.cores = 8)

# Comparison of AIC

sum(result[1, ] < result[2, ]) / N.sim
sum(result[1, ] < result[3, ]) / N.sim
sum(result[1, ] < result[4, ]) / N.sim

# Coefficients in the path analysis

## y ~ x
summary(result[5, ])

## y ~ z1
summary(result[6, ])