Created
May 1, 2020 09:51
-
-
Save vanAmsterdam/40f6213b5827f3b1ce9d2772f1c7f53f to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #' 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 | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment