latentconfoundermodel.R

#' Latent confounder model

# dag: 

#' W1 <- U -> W2
#' tx -> y
#' U -> y
#' U -> tx

library(lavaan)

## simulate data
set.seed(1234)
nsim = 1000

s_W1    = 0.1
s_W2    = 0.2
s_tx   = 0.2
s_y    = 0.2
b_tx_y = 1
b_U_W1 = 0.5
b_U_W2 = 0.5
b_U_tx = .75
b_U_y  = .75

U  = rnorm(nsim)
W1 = U * b_U_W1 +  rnorm(nsim, sd=s_W1)
W2 = U * b_U_W2 +  rnorm(nsim, sd=s_W2)
tx = U * b_U_tx + rnorm(nsim, sd=s_tx)
### create potential outcomes and observed outcome
y0 = U * b_U_y  + rnorm(nsim, sd=s_y)
y1 = y0 + b_tx_y
y  = y0 + b_tx_y * tx
dfall       = data.frame(U, W1, W2, tx, y, y0, y1)
dfall$train = rep(c(T,F), each=nsim / 2)
dftrain     = dfall[dfall$train==T,setdiff(colnames(dfall), 'U')]
dftest      = dfall[dfall$train==F,setdiff(colnames(dfall), 'U')]
dftestdotx  = dftest
dftestdotx$tx = rnorm(nsim/2, mean(mean(tx), sd=sd(tx)))
dftestdotx$tx = rbinom(nsim/2, size=1, prob=0.5)
dftestdotx$y  = ifelse(dftestdotx$tx, dftestdotx$y1, dftestdotx$y0)

## make groundtruth model
fittrue = lm(y~U+tx)
summary(fittrue)

## flat model
fitflat = lm(y~tx+W1+W2, data=dftrain)
summary(fitflat)
dftest$yhatflat     = predict(fitflat, newdata=dftest)
dftestdotx$yhatflat = predict(fitflat, newdata=dftestdotx)

## SEM model with lavaan
m1 = '
  # regressions
  y ~ b_tx_y * tx
  
  # latents
  U =~ NA*W1 + b_U_W1*W1 + b_U_W2*W2 + b_U_tx * tx + b_U_y * y # using NA before W makes lavaan regress this factor instead of fixing it to 1
  
  # constraints
  U ~~ 1 * U # fixes variance to 1
  b_U_tx > 0
  b_U_W1 > 0
  b_U_W2 > 0
  b_U_y  > 0
'
fs1 = sem(m1, data=dftrain)
coef(fs1)

# model implied covariance matrix of observed variables
cv <- fitted(fs1)$cov
cv