Set up our causal model.

n <- 1e2
U <- rnorm(n)
T <- 4 * U + rnorm(n)
Y <- 8 * T -5 * U + rnorm(n)
O <- as.numeric(cut(U, 5))

dat <- data.frame(cbind(U, T, Y, O))

Look at our Naive estimate

lm(Y ~ T + O, dat)
## 
## Call:
## lm(formula = Y ~ T + O, data = dat)
## 
## Coefficients:
## (Intercept)            T            O  
##       7.581        7.375       -2.581

Look at our bayesian bootstrapping method

library(ordinalconfounder)
print(estimate(dat))
## Running MCMC with 4 parallel chains...
## 
## Chain 2 finished in 4.8 seconds.
## Chain 3 finished in 5.1 seconds.
## Chain 4 finished in 5.1 seconds.
## Chain 1 finished in 5.4 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 5.1 seconds.
## Total execution time: 5.5 seconds.
## 
## Running MCMC with 4 parallel chains...
## 
## Chain 1 finished in 6.8 seconds.
## Chain 2 finished in 7.1 seconds.
## Chain 4 finished in 7.1 seconds.
## Chain 3 finished in 7.2 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 7.1 seconds.
## Total execution time: 7.3 seconds.
## 
##        T 
## 7.896396