Elastic Net Cp function

enet_cp.r

#sigma2, df, and Cp may not be reliable in default elasticnet function enet
#provide replacements
#S2L is the estimate for sigma2 based on maximal model (full_ols_model from lm())
S2L<- sum(residuals(full_ols_model)^2)/full_ols_model$df.residual

Cp_enet<-function(enet_obj,y,X,S2L){
  #returns the Mallows Cp statistics for a given elastic net object
  betamat<-enet_obj$beta.pure
  #create a matrix where each column is prediction at one of the df levels
  yhats<-X%*%t(betamat) + mean(y)
  resids<-yhats-y
  ESS<-colSums(resids^2)
  lambda<-enet_obj$lambda
  #indices of which columns are in active set
  actives<-unlist(enet_obj$actions) 
  k<-ncol(X)
  df<-c(1,rep(NA,k))
  #df is determined by trace of hat matrix
  for(i in 1:k){
    Xactive<-X[,actives[1:i]]
    H = Xactive%*%solve(t(Xactive)%*%Xactive+diag(lambda,i),t(Xactive))
    df[i+1] = 1+sum(diag(H)) #compute trace
  }
  n<-length(y)
  Cp<-ESS/S2L-n+2*df
  return(Cp) #a vector
}