Today's lab is about implementing a Logistic Regression with a Ridge Regression penalty using Newton's Method.

You may use these notes provided by the TA for the previous semester for implementing your code. logistic regression

The data set for the lab was taken from “Regression With Graphics” by Hamilton. A town was considering closing a school due to toxic material on the grounds, and while the debate was on going, researchers collected data about whether individuals in town favored the closing along with other covariates such as years they lived in town, whether they have children, etc.

We will focus in on the “favor” variable as the response, which is 0 if they did not favor closing the school and 1 if they did.

And for the covariate, we will consider the years they lived in the town which is an integer on 1 to 81.

We seek to model the probability that an individual will favor closing the school given the years they lived in the town.

Exercises

Implementing Exercise 1 is your goal for the lab, but first load the data. Get the data from here and save it in your working directory(even though this file has an R file extension it is actually a tab delimited table). Load the data with the following R commands:

   data = read.table("close.school.for.r",header=TRUE)
   #The bernoulli response "favor"
   z = data[,1]
   #A column of intercepts and the years lived in town.
   X = cbind(1,data[,2])
   X = as.matrix(X)
   colnames(X) = c("intercept", "years")

1. Write functions to calculate penalized maximum likelihood estimates for a logistic regression using Newton's Method. Create estimates of parameters for intercepts and years with favor closing (z) as the response.

gradient = function(z, X, beta, lambda){
p = pi(X, beta)
g = -t(X) %*% (z - p) + lambda*beta
#print(g)
return(g)
 }
# hessian 
# L0 is lambda*I_0 in the notes.
hessian = function(X, beta, lambda){
I = diag(ncol(X))
p = pi(X, beta)
W = diag(p*(1-p))
#print(lambda*I)
H = t(X) %*% W %*% X + lambda*I
#print(H)
return(H)
}

2. Verify that your MLEs when lambda = 0 are the same estimates as those produced by Rs built in function to do a logistic regression. Use the following code in R to get those estimates:

fit = glm(z ~ X-1, family="binomial")
coef(fit)

3. Plot contours of the penalized objective function around values of the fit parameters with the following code:

pi = function(X, beta){
return(as.vector(1/(1+exp(- X %*% beta))))
}
objective = function(z, X, beta, lambda){
p = pi(X,beta)
return(-sum(z*log(p) + (1-z)*log(1-p)) + lambda/2*sum(beta^2))
}
intercepts = seq(-2,2,,100)
yearcoefs = seq(-.25,.25,,100)
grid = expand.grid(intercepts,yearcoefs)
grid = as.matrix(grid)
obj = apply(grid,1, function(b) objective(z,X,b,lambda))
obj = matrix(obj,100,100)
zlim = quantile(obj,c(0,.65))  #verify is finite
contour(intercepts,yearcoefs,obj, zlim=zlim,xlab="intercept", ylab="year coef")

4. Repeat the same using Gradient descent and conjugate gradient method.

lab13.r