Logistic Regression: STAT 963

Step 1: Load and clean up data

# Load the data
dat <- ElemStatLearn::SAheart
vo <- c('sbp','tobacco','ldl','famhist','obesity','alcohol','age')
X <- dat[,vo]
# Re-code famhist
X <- cbind(intercept=1,famhist=model.matrix(~-1+famhist,data=dat)[,2],X[,-which(colnames(X)=='famhist')])
X <- as.matrix(X)[,c('intercept',vo)]
y <- dat$chd
# Define the p(x_i;beta) function
p <- function(xi,beta) { 1/(1+exp(-(sum(xi*beta)))) }
head(cbind(y,X))

  y intercept sbp tobacco  ldl famhist obesity alcohol age
1 1         1 160   12.00 5.73       1   25.30   97.20  52
2 1         1 144    0.01 4.41       0   28.87    2.06  63
3 0         1 118    0.08 3.48       1   29.14    3.81  46
4 1         1 170    7.50 6.41       1   31.99   24.26  58
5 1         1 134   13.60 3.50       1   25.99   57.34  49
6 0         1 132    6.20 6.47       1   30.77   14.14  45

Step 2: Estimate coefficients with Newton-Raphson

# Set up optimization parameters
beta.old <- rep(0,ncol(X))
tol <- 10e-3
dist <- 1
i <- 0
# Begin the loop!
while(dist > tol) {
  i <- i + 1
  # Get the p-vector
  p.vec <- apply(X,1,p,beta=beta.old)
  # Get the W weight matrix
  W <- diag(p.vec*(1-p.vec))
  # Get the Score vector (i.e. gradient)
  dldbeta <- t(X) %*% (y-p.vec)
  # Get the Hessian
  dl2dbeta2 <- -t(X) %*% W %*% X
  # Update beta
  beta.new <- beta.old - solve(dl2dbeta2) %*% dldbeta
  # Calculate distance
  dist <- sqrt(sum((beta.new-beta.old)^2))
  # Reset
  beta.old <- beta.new
}
paste('Convergence achieved after',i,'steps',sep=' ')

[1] "Convergence achieved after 4 steps"

round(beta.new,3)

            [,1]
intercept -4.130
sbp        0.006
tobacco    0.080
ldl        0.185
famhist    0.939
obesity   -0.035
alcohol    0.001
age        0.043

Step 3: Estimate coefficients with IRLS algorithm

# Set up optimization parameters
beta.old <- rep(0,ncol(X))
tol <- 10e-3
dist <- 1
i <- 0
# Begin the loop!
while(dist > tol) {
  i <- i + 1
  # Get the p-vector
  p.vec <- apply(X,1,p,beta=beta.old)
  # Get the W weight matrix
  W <- diag(p.vec*(1-p.vec))
  # Get W inverse
  Winv <- diag(1/(p.vec*(1-p.vec)))
  # Get z
  z <- (X %*% beta.old) + Winv %*% (y-p.vec)
  # Update beta
  beta.new <- solve(t(X) %*% W %*% X) %*% t(X) %*% W %*% z
  # Calculate distance
  dist <- sqrt(sum((beta.new-beta.old)^2))
  # Reset
  beta.old <- beta.new
}
paste('Convergence achieved after',i,'steps',sep=' ')

[1] "Convergence achieved after 4 steps"

round(beta.new,3)

            [,1]
intercept -4.130
sbp        0.006
tobacco    0.080
ldl        0.185
famhist    0.939
obesity   -0.035
alcohol    0.001
age        0.043