Bayesian implementation of instrumental variable regression (using rstan)

bayes_instrument.R

# install rethinking library with instructions here:
# https://github.com/rmcelreath/statrethinking_winter2019

library(rethinking)
library(dagitty)
library(shape)

# the structural model
# W: wages
# Q: quarter of year of birth
# E: education
# U: unobserved confound

the_dag <- dagitty("dag{
    Q -> E -> W
    E <- U -> W
}")
coordinates(the_dag) <- list(x=c(Q=0,E=1,W=2,U=1.5),y=c(U=0,Q=1,E=1,W=1))
drawdag(the_dag)

# find instruments
instrumentalVariables( the_dag , exposure="E" , outcome="W" )

# simulate
# bEW: effect of education on wages

bEW <- 0

set.seed(73)
N <- 500
U_sim <- rnorm( N )
Q_sim <- sample( 1:4 , size=N , replace=TRUE )
E_sim <- rnorm( N , U_sim + Q_sim )
W_sim <- rnorm( N , U_sim + bEW*E_sim )
dat_sim <- list( 
    W=standardize(W_sim) , 
    E=standardize(E_sim) , 
    Q=standardize(Q_sim) )

# standard regression of W on E (is biased by U)
m14.4 <- ulam(
    alist(
        W ~ dnorm( mu , sigma ),
        mu <- aW + bEW*E,
        aW ~ dnorm( 0 , 0.2 ),
        bEW ~ dnorm( 0 , 0.5 ),
        sigma ~ dexp( 1 )
    ) , data=dat_sim , chains=4 , cores=4 )
precis( m14.4 )

# instrumental variable regression (unbiased)
# note covariance between W and E measured by Rho[1,2]

m14.5 <- ulam(
    alist(
        c(W,E) ~ multi_normal( c(muW,muE) , Rho , Sigma ),
        muW <- aW + bEW*E,
        muE <- aE + bQE*Q,
        c(aW,aE) ~ normal( 0 , 0.2 ),
        c(bEW,bQE) ~ normal( 0 , 0.5 ),
        Rho ~ lkj_corr( 2 ),
        Sigma ~ exponential( 1 )
    ), data=dat_sim , chains=4 , cores=4 )
precis( m14.5 , depth=3 )