19 Ridge Regression
THIS CHAPTER IS UNDER CONSTRUCTION!!!
We should provide an example in Stan.
19.1 Introduction
# Settings
library(R2OpenBUGS)
<- "C:/Program Files/OpenBUGS323/OpenBugs.exe" # location of OpenBUGS
bugslocation <- file.path(getwd(), "BUGS") # Bugs working directory
bugsworkingdir
#-------------------------------------------------------------------------------
# Simulate fake data
#-------------------------------------------------------------------------------
library(MASS)
<- 50 # sample size
n <- 1.2
b0 <- rnorm(5, 0, 2)
b <- matrix(c(10,3,3,2,1,
Sigma 3,2,3,2,1,
3,3,5,3,2,
2,2,3,10,3,
1,1,2,3,15),5,5)
Sigma<- mvrnorm(n = n, rep(0, 5), Sigma)
x
<- rnorm(n, 0, sd=3) # residuals
simresid
<- x
x.z for(i in 1:ncol(x)) x.z[,i] <- (x[,i]-mean(x[,i]))/sd(x[,i])
<- b0 + x.z%*%b + simresid # calculate y, i.e. the data
y
#-------------------------------------------------------------------------------
# Function to generate initial values
#-------------------------------------------------------------------------------
<- function() {
inits list(b0=runif(1, -2, 2),
b=runif(5, -2, 2),
sigma=runif(1, 0.1, 2))
}
#-------------------------------------------------------------------------------
# Run OpenBUGS
#-------------------------------------------------------------------------------
<- c("b0", "b", "sigma")
parameters
<- c(1, 2, 10, 25, 50, 100, 500, 1000, 10000)
lambda <- matrix(ncol=length(lambda), nrow=length(b))
bs <- matrix(ncol=length(lambda), nrow=length(b))
bse
for(j in 1:length(lambda)){
<- list(y=as.numeric(y), x=x, n=n, mb=rep(0, 5), lambda=lambda[j])
datax
<- bugs(datax, inits, parameters, model.file="ridge_regression.txt",
fit n.thin=1, n.chains=2, n.burnin=5000, n.iter=10000,
debug=FALSE, OpenBUGS.pgm = bugslocation,
working.directory=bugsworkingdir)
<- fit$mean$b
bs[,j] <- fit$sd$b
bse[,j]
}
range(bs)
plot(1:length(lambda), seq(-2, 1, length=length(lambda)), type="n")
<- rainbow(length(b))
colkey for(j in 1:nrow(bs)){
lines(1:length(lambda), bs[j,], col=colkey[j], lwd=2)
lines(1:length(lambda), bs[j,]-2*bse[j,], col=colkey[j], lty=3)
lines(1:length(lambda), bs[j,]+2*bse[j,], col=colkey[j], lty=3)
}abline(h=0)
round(fit$summary,2)
#-------------------------------------------------------------------------------
# Run WinBUGS
#-------------------------------------------------------------------------------
library(R2WinBUGS)
<- "C:/Users/fk/WinBUGS14" #
bugsdir <- bugs(datax, inits= inits, parameters,
mod model.file="normlinreg.txt", n.chains=2, n.iter=1000,
n.burnin=500, n.thin=1, debug=TRUE,
bugs.directory=bugsdir, program="WinBUGS", working.directory=bugsworkingdir)
#-------------------------------------------------------------------------------
# Test convergence and make inference
#-------------------------------------------------------------------------------
library(blmeco)
# Make Figure 12.2
par(mfrow=c(3,1))
historyplot(fit, "beta0")
historyplot(fit, "beta1")
historyplot(fit, "sigmaRes")
# Parameter estimates
print(fit$summary, 3)
# Make predictions for covariate values between 10 and 30
<- data.frame(x=seq(10, 30, length=100))
newdat <- model.matrix(~x, data=newdat)
Xmat <- matrix(ncol=fit$n.sim, nrow=nrow(newdat))
predmat for(i in 1:fit$n.sim) predmat[,i] <- Xmat%*%c(fit$sims.list$beta0[i], fit$sims.list$beta1[i])
$lower.bugs <- apply(predmat, 1, quantile, prob=0.025)
newdat$upper.bugs <- apply(predmat, 1, quantile, prob=0.975)
newdat
plot(y~x, pch=16, las=1, cex.lab=1.4, cex.axis=1.2, type="n", main="")
polygon(c(newdat$x, rev(newdat$x)), c(newdat$lower.bugs, rev(newdat$upper.bugs)), col=grey(0.7), border=NA)
abline(c(fit$mean$beta0, fit$mean$beta1), lwd=2)
box()
points(x,y)