ridge_regression.R

#### Initialization
# set params
stan.chains=4
stan.iter=10000
stan.thin=5
max_treedepth=20
adapt_delta=0.9

# load deps
library(glmnet)
library(rstan)
library(shinystan)

# global pars
options(mc.cores=4)
rstan_options(auto_write=TRUE)

#### Preliminary processing
# prepare data
data(BinomialExample)
x = model.matrix(y ~ (.), data=data.frame(cbind(y,x)))[,-1]

#### Main processing
model_code="
	data {
		int p; // number covariates
		int n; // number covariates
		matrix[n,p] x; // number covariates
		int y[n]; // response
	}
	
	transformed data {
		// z-score data
		vector[p] x_mean;
		vector[p] x_sd;
		matrix[n,p] x_std;
		
		for (i in 1:p) {
			x_mean[i] <- mean(col(x, i));
			x_sd[i] <- sd(col(x, i));
			for (j in 1:n)
				x_std[j,i] <- (x[j,i] - x_mean[i])/ x_sd[i];
		}
	}
	
	parameters {
		real beta0;
		vector[p] beta;
		real<lower=0> phi;
	}
	
	transformed parameters {
		vector[n] mu;
		mu <- beta0 + x_std * beta;
	}
	
	model {
		// priors
		beta0 ~ normal(0, 10);
		phi ~ gamma(0.1, 0.1);
		beta ~ normal(0, sqrt(phi));
		
		// likelihood
		y ~ bernoulli_logit(mu);
	}

"


model=stan(
	model_code=model_code,
	data=list(
		x=x,
		y=y,
		p=ncol(x),
		n=nrow(x)
	),
	chains=stan.chains,
	iter=stan.iter,
	thin=stan.thin,
	control=list(
		max_treedepth=max_treedepth,
		adapt_delta=adapt_delta
	)
)

#### Post analysis
launch_shinystan(model)