2.11-ridge_regression.Rmd

# Ridge Regression {#ridge_regression}

THIS CHAPTER IS UNDER CONSTRUCTION!!! 

> We should provide an example in Stan.

## Introduction


```{r, eval=FALSE}
# Settings 
library(R2OpenBUGS)
bugslocation <- "C:/Program Files/OpenBUGS323/OpenBugs.exe"      # location of OpenBUGS
bugsworkingdir <- file.path(getwd(), "BUGS")                     # Bugs working directory

#-------------------------------------------------------------------------------
# Simulate fake data 
#-------------------------------------------------------------------------------
library(MASS)
n <- 50                                      # sample size
b0 <- 1.2
b <- rnorm(5, 0, 2)
Sigma <- matrix(c(10,3,3,2,1,
                  3,2,3,2,1,
                  3,3,5,3,2,
                  2,2,3,10,3,
                  1,1,2,3,15),5,5)
Sigma
x <- mvrnorm(n = n, rep(0, 5), Sigma)

simresid <- rnorm(n, 0, sd=3)            # residuals

x.z <- x
for(i in 1:ncol(x)) x.z[,i] <- (x[,i]-mean(x[,i]))/sd(x[,i]) 

y <- b0 + x.z%*%b + simresid                    # calculate y, i.e. the data



#-------------------------------------------------------------------------------
# Function to generate initial values 
#-------------------------------------------------------------------------------
inits <- function() {
  list(b0=runif(1, -2, 2),
       b=runif(5, -2, 2),
       sigma=runif(1, 0.1, 2))
}


#-------------------------------------------------------------------------------
# Run OpenBUGS
#-------------------------------------------------------------------------------
parameters <- c("b0", "b", "sigma")

lambda <- c(1, 2, 10, 25, 50, 100, 500, 1000, 10000)
bs <- matrix(ncol=length(lambda), nrow=length(b))
bse <- matrix(ncol=length(lambda), nrow=length(b))

for(j in 1:length(lambda)){
datax <- list(y=as.numeric(y), x=x, n=n, mb=rep(0, 5), lambda=lambda[j])


fit <- bugs(datax, inits, parameters, model.file="ridge_regression.txt",
            n.thin=1, n.chains=2, n.burnin=5000, n.iter=10000,
            debug=FALSE, OpenBUGS.pgm = bugslocation, 
            working.directory=bugsworkingdir)

bs[,j] <- fit$mean$b
bse[,j] <- fit$sd$b

}

range(bs)
plot(1:length(lambda), seq(-2, 1, length=length(lambda)), type="n")
colkey <- rainbow(length(b))
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)
bugsdir <- "C:/Users/fk/WinBUGS14"         # 
mod <- bugs(datax, inits= inits, parameters,
            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
newdat <- data.frame(x=seq(10, 30, length=100))
Xmat <- model.matrix(~x, data=newdat)
predmat <- matrix(ncol=fit$n.sim, nrow=nrow(newdat))
for(i in 1:fit$n.sim) predmat[,i] <- Xmat%*%c(fit$sims.list$beta0[i], fit$sims.list$beta1[i])
newdat$lower.bugs <- apply(predmat, 1, quantile, prob=0.025)
newdat$upper.bugs <- apply(predmat, 1, quantile, prob=0.975)

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)
```