simulaton of instrumental variable regression with a count outcome

Count_IV.R

library(rethinking)

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 + 0*E_sim )

Rsize <- 2
R_sim <- rbinom( N , size=Rsize , inv_logit(W_sim) )

dat_sim <- list(
    W=standardize(W_sim) ,
    E=standardize(E_sim) ,
    Q=standardize(Q_sim) ,
    R=R_sim,
    Rsize=Rsize )

m14.6 <- 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 )

mc <- "
data{
    array[500] int R;
     vector[500] W;
     vector[500] E;
     vector[500] Q;
     int Rsize;
}
parameters{
     real aE;
     real aW;
     real bQE;
     real bEW;
     corr_matrix[2] Rho;
     vector<lower=0>[2] Sigma;
     vector[500] logit_p;
}
model{
     vector[500] muW;
     vector[500] muE;
    Sigma ~ exponential( 1 );
    Rho ~ lkj_corr( 2 );
    bEW ~ normal( 0 , 0.5 );
    bQE ~ normal( 0 , 0.5 );
    aW ~ normal( 0 , 0.2 );
    aE ~ normal( 0 , 0.2 );
    for ( i in 1:500 ) {
        muE[i] = aE + bQE * Q[i];
    }
    for ( i in 1:500 ) {
        muW[i] = aW + bEW * E[i];
    }
    {
    array[500] vector[2] YY;
    array[500] vector[2] MU;
    for ( j in 1:500 ) MU[j] = [ muW[j] , muE[j] ]';
    for ( j in 1:500 ) YY[j] = [ logit_p[j] , E[j] ]';
    YY ~ multi_normal( MU , quad_form_diag(Rho , Sigma) );
    }
    R ~ binomial_logit(Rsize,logit_p);
}
"

m <- cstan( model_code=mc , data=dat_sim , chains=1 , rstan_out=FALSE )

precis(m)

precis(m14.6)