Navigation Menu

Skip to content

Instantly share code, notes, and snippets.

@rmcelreath
Created June 29, 2021 10:57
Show Gist options
  • Star 1 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save rmcelreath/9c687b2cf504c0c2910673a056ffde45 to your computer and use it in GitHub Desktop.
Save rmcelreath/9c687b2cf504c0c2910673a056ffde45 to your computer and use it in GitHub Desktop.
Code eamples for RFDT
# more exotic example - no instrument (B1 -> D), but have two measures of U
set.seed(1908)
N <- 200 # number of pairs
U <- rnorm(N,0,1) # simulate confound
V <- rnorm(N,U,1)
W <- rnorm(N,U,1)
# birth order and family sizes
B1 <- rbinom(N,size=1,prob=0.5) # 50% first borns
M <- rnorm( N , 2*B1 + U )
B2 <- rbinom(N,size=1,prob=0.5)
D <- rnorm( N , 2*B2 + 0.5*B1 + U + 0*M )
# confounded regression
precis( lm( D ~ M + B1 + B2 + V + W ) )
# full-luxury bayesian inference
dat2 <- list(N=N,M=M,D=D,B1=B1,B2=B2,V=V,W=W)
flbi2 <- ulam(
alist(
M ~ normal( muM , sigmaM ),
muM <- a1 + b*B1 + k*U[i],
D ~ normal( muD , sigmaD ),
muD <- a2 + b*B2 + d*B1 + m*M + k*U[i],
W ~ normal( muW , sigmaW ),
muW <- a3 + w*U[i],
V ~ normal( muV , sigmaV ),
muV <- a4 + v*U[i],
vector[N]:U ~ normal(0,1),
c(a1,a2,a3,a4,b,d,m) ~ normal( 0 , 0.5 ),
c(k,w,v) ~ exponential( 1 ),
c(sigmaM,sigmaD,sigmaW,sigmaV) ~ exponential( 1 )
), data=dat2 , chains=4 , cores=4 , iter=2000 , cmdstan=TRUE )
precis(flbi2)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment