Proxy problem: https://twitter.com/y2silence/status/1251151157264674816

proxy.R

library(rstan)

stan_code <- "
 data {
  int<lower=0> N;
  vector[N] y;
  vector[N] a;
  vector[N] x1;
  vector[N] x2;
}

parameters {
  vector[N] u;
  real b_x1; // note potential multimodality due to the switching signs of b_x1 and b_x2
  real b_x2; // could use ordered[2] for these...
  real b_a;
  real b_ya;
  real b_yu;
  real a_x1;
  real a_x2;
  real a_a;
  real a_y;
  real<lower=0> sigma_x1;
  real<lower=0> sigma_x2;
  real<lower=0> sigma_a;
  real<lower=0> sigma_y;
}

model {
  a_x1 ~ normal(0, 1);
  a_x2 ~ normal(0, 1);
  a_a ~ normal(0, 1);
  a_y ~ normal(0, 1);
  b_x1  ~ normal(0, 1);
  b_x2  ~ normal(0, 1);
  b_a  ~ normal(0, 1);
  b_ya  ~ normal(0, 1);
  b_yu  ~ normal(0, 1);
  u ~ normal(0, 1); // prior directly on U
  x1 ~ normal(a_x1 + b_x1 * u, sigma_x1);
  x2 ~ normal(a_x2 + b_x2 * u, sigma_x2);
  a ~ normal(a_a + b_a * u, sigma_a);
  y ~ normal(a_y + b_ya * a + b_yu * u, sigma_y);
}
"

set.seed(1)
n <- 1000

u <- rnorm(n, 0, 1)
x1 <- rnorm(n, -u, 1)
x2 <- rnorm(n, 2*u, 1)
a <- rnorm(n, u, 1)
y <- rnorm(n, u + a, 1)

fit <- stan(model_code = stan_code, data = tibble::lst(a,y,x1,x2,N=n), pars = "u", include = FALSE,
                chains = 4, cores = 4, iter = 5000)

print(fit)