302_ridge.Rmd

# Ridge Regression - L2 Regularization

<!-- Discussion of the model goes here -->

```{r}

ridge <- function(X, Y, rm_na=T, standardise=F, lambda, method="matrix") {
  
  if (!is.matrix(X)) {
    X <- as.matrix(X)
  }
  
  if (rm_na) {
    index <- complete.cases(X) & complete.cases(Y)
    X <- X[index,]
    Y <- Y[index]
  }
  
  if (standardise) {
    X <- apply(X, 2, function(x) (x - mean(x))/sd(x))
  }
  
  # if (intercept) {
  #   X <- cbind(rep(1,nrow(X)), X)
  #   colnames(X)[1] <- "(Intercept)"
  # }
  
  k <- ncol(X)
  
  if (method=="matrix") {
    b <- solve(crossprod(X) + lambda*diag(k)) %*% crossprod(X, Y)
  } else if (method=="MLE") {
    start <- numeric(k)
    ridgeMin <- function(b) {
      crossprod(Y - X%*%b) + lambda * crossprod(b)
    }
    out <- optim(start, ridgeMin, method="BFGS")
    b <- out$par
    names(b) <- colnames(X)
  }
  
  return(b)
  
}


```

Testing the function:

```{r}

library(glmnet)

X <- as.matrix(mtcars[, -1])
X_standard <- apply(X, 2, function(x) (x-mean(x)) / sd(x))
Y <- mtcars[[1]]

data.frame(
  glmnet = glmnet(X, Y, alpha=0, lambda=0.5, intercept = F, standardize=F)$beta[,1],
  Ridge_Matrix = ridge(X, Y, lambda=0.5),
  Ridge_MLE = ridge(X, Y, lambda=0.5, method="MLE")
)

data.frame(
  glmnet = glmnet(X_standard, Y, alpha=0, lambda=0.5, intercept = F, standardize=F)$beta[,1],
  Ridge_Matrix = ridge(X_standard, Y, lambda=0.5),
  Ridge_MLE = ridge(X_standard, Y, lambda=0.5, method="MLE")
)

```