# Regression example, taken from “A First Course in Bayesian Statistical
# Methods” by Peter D. Hoff
# Arizona diabetes data
rm(list=ls()); set.seed(20210325); S <- 5000; library(MASS); # 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); round(apply(X, 2, mean)); apply(X, 2, sd); X <- cbind(1, X); colnames(X) <- c("ones", names); rm(names); # 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) # Prelim calculations n <- dim(X)[1]; k <- dim(X)[2]; beta.hat <- solve( t(X) %*% X, t(X) %*% y )[,1]; beta.hat; XTX.inv <- solve( t(X) %*% X ); XTX.inv; y.hat <- (X %*% beta.hat)[,1]; # y.hat; s2 <- sum((y - y.hat)^2)/(n-k); s2; rm(y.hat); # Check the matrix calculations versus lm() results m0 <- lm(y ~ 0 + X); coef(m0) # should be beta.hat sigma(m0)^2 # should be s2 vcov(m0) / sigma(m0)^2 # Should be XTX.inv rm(m0) # Posterior sampling sigma2.chain <- rep(NA, S); beta.chain <- matrix(rep(NA, S*k), S, k); for(s in 1:S) { sigma2 <- (n-k)*s2 / rchisq(1, df=n-k) beta <- mvrnorm(1, mu=beta.hat, Sigma=sigma2*XTX.inv) sigma2.chain[s] <- sigma2 beta.chain[s,] <- beta } sigma.chain <- sqrt(sigma2.chain) probs <- c(.025, .25, .5, .75, .975) Results <- apply(beta.chain, 2, quantile, probs=probs) Results <- t(Results); Results <- rbind(Results, quantile(sigma.chain, probs=probs)) rownames(Results) <- c(colnames(X), "sigma") round(Results, 2) op <- par(mfrow=c(2,4)) for(j in 1:7) { hist(beta.chain[,j], freq=F, xlab="", ylab="", breaks=30, main=paste("Posterior dist of coefficient of", rownames(Results)[j])); lines(density(beta.chain[,j]), col="red") } hist(sigma.chain, freq=F, xlab="", ylab="", breaks=30, main="Posterior dist for error SD"); lines(density(sigma.chain), col="red") par(op) # Predictions?!?!? # Okay, let's suppose we wish to predict the response (glucose) for # an individual whose x-values are given by x.tilde <- X[21,]; x.tilde; # Easy enough y.tilde <- rnorm(S, mean=(beta.chain%*%x.tilde)[,1], sd=sigma.chain) hist(y.tilde, freq=F, breaks=30, xlab="", ylab="", main="Posterior predictive dist for y.tilde") lines(density(y.tilde), col="red") quantile(y.tilde, probs=probs)