mid

rm(list=ls())
library(glmnet); library(ncvreg); library(MASS)

## Loading required package: Matrix

## Loading required package: foreach

## Loaded glmnet 2.0-18

e.fun = function(y.vec,x.mat,g.vec,n.mod=5){
  b.mat = matrix(0,ncol(x.mat),n.mod)
  #xy.df = data.frame(y.vec,x.mat)
  #n.mod = lm(y.vec~1,data=xy.df); f.mod = lm(y.vec~.,data=xy.df)
  ret = NULL
  for(g in g.vec){
    fit = cv.ncvreg(x.mat,y.vec,gamma=g)
    opt = which.min(fit$cve)
    ret = cbind(ret,c(fit$lambda[opt],fit$cve[opt]))
  }
  opt = which.min(ret[2,])
  opt.parm = c(ret[1,opt],g.vec[opt])
  b.mat[,1] = coef(ncvreg(x.mat,y.vec,lambda=opt.parm[1],gamma=opt.parm[2]))[-1] # SCAD1
  b.mat[,2] = coef(cv.glmnet(x.mat,y.vec))[-1] #Lasso
  b.mat[,3] = coef(cv.glmnet(x.mat,y.vec,alpha=0,))[-1] #Ridge
  b.mat[,4] = coef(lm(y.vec~x.mat))[-1] #OLS
  b.mat[tb.vec!=0,5] = coef(lm(y.vec~x.mat[,tb.vec!=0]))[-1] #Oracle
  
  return(b.mat)
}

m.fun = function(b.mat,tb.vec){
  sen = colMeans(b.mat[tb.vec!=0,]!=0)
  spc = ifelse(sum(tb.vec==0)==0,0,colMeans(b.mat[tb.vec==0,]==0))
  acc = (sen+spc)==2
  mse = colSums((b.mat-tb.vec)^2)
  return(cbind(sen,spc,acc,mse))
}

Including Plots

You can also embed plots, for example: