# Regression example, taken from “A First Course in Bayesian Statistical
# Methods” by Peter D. Hoff
# Arizona diabetes data
rm(list=ls()); set.seed(20210325);
# Read in data
file <- "http://www2.stat.duke.edu/~pdh10/FCBS/Exercises/azdiabetes.dat"
Data <- read.table(file=file, header=T); rm(file);
dim(Data); names(Data);
# Data pre-processing
y <- Data$glu;
X <- as.matrix(Data[,c(1,3,4,5,6,7)]); rm(Data);
names <- colnames(X);
m <- apply(X, 2, mean); s <- apply(X, 2, sd);
X <- sweep(X, 2, m, "-"); rm(m);
X <- sweep(X, 2, s, "/"); rm(s);
apply(X, 2, mean); apply(X, 2, sd);
X <- cbind(1, X); colnames(X) <- c("ones", names);
n <- dim(X)[1]; k <- dim(X)[2];
# We assume the standard noninformative prior on the variance,
# p(beta, sigma2) = const*(sigma2)^(-1)
# This gives the marginal posterior
# sigma2 | y ~ Inv-Chisq(n-k, s2)
# where s2 = 1/(n-k) * (y - X %*% beta.hat)' (y - X %*% beta.hat)
# and
# beta | sigma2, y ~ Normal with mean vector and covariance matrix:
# (X'X)^(-1) %*% X' y & sigma2 * (X'X)^(-1)
library(rstan)
stan_model <- "
data{
int
int
vector [n] y;
matrix[n,k] X;
# vector[k] beta0;
# matrix[k,k] Omega;
}
parameters{
vector[k] beta;
# real
real log_sigma;
}
transformed parameters{
vector[n] mu;
real
mu = X * beta;
sigma <- exp(log_sigma); # sqrt(sigma2);
}
model{
# beta ~ multi_normal(beta0, sigma2*Omega) ;
# sigma2 ~ inv_gamma(.001, .001);
y ~ normal(mu, sigma) ;
}
"
beta0 <- rep(0, k)
data <- list(n=n, k=k, y=y, X=X) # , beta0=beta0, Omega=Omega)
fit <- stan(model_code=stan_model,
data=data, iter=2000, chains=5)
Results <- summary(fit)$summary[c(1:7, 541), ]
Results
round(Results[, 4:8], 2)
sims <- extract(fit)