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)
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))
}
set.seed(201972026)
n.vec = c(50,100,200,400,800,1600)
g.vec = 2*(1:10)+0.1
p = c(10,100)
n.mod = 5
Sensitivity = Specificity = Accuracry = MSE = NULL
iters=100
for(m in p){
tbvec = list(seq(1,m),c(seq(1,5),rep(0,m-5)))
SN=0
for (tb.vec in tbvec){
for (i in 1:length(tb.vec)){tb.vec[i]=ifelse(tb.vec[i]!=0,1/tb.vec[i],tb.vec[i])}
SN = SN + 1
bn = length(tb.vec)
mu = rep(0, bn)
for(n in 1:length(n.vec)){
sen.list = spc.list = acc.list = mse.list = list()
for(iter in 1:iters){
corr.mat = outer(1:bn, 1:bn, function(x,y) {0.5^abs(x-y)})
x.mat = mvrnorm(n.vec[n], mu = mu, Sigma = corr.mat)
y.vec = drop(x.mat%*%drop(tb.vec)+rnorm(n.vec[n],0,1))
b.mat = e.fun(y.vec,x.mat,g.vec)
sen.list[[iter]] = m.fun(b.mat,tb.vec)[,1]
spc.list[[iter]] = m.fun(b.mat,tb.vec)[,2]
acc.list[[iter]] = m.fun(b.mat,tb.vec)[,3]
mse.list[[iter]] = m.fun(b.mat,tb.vec)[,4]
}
Sensitivity = rbind(Sensitivity,cbind(SN,m,n.vec[n],t(Reduce("+",sen.list)/iters)))
Specificity = rbind(Specificity,cbind(SN,m,n.vec[n],t(Reduce("+",spc.list)/iters)))
Accuracry = rbind(Accuracry,cbind(SN,m,n.vec[n],t(Reduce("+",acc.list)/iters)))
MSE = rbind(MSE,cbind(SN,m,n.vec[n],t(Reduce("+",mse.list)/iters)))
}
}
}
colnames(Sensitivity) = c('beta.vec', 'p', 'N', 'SCAD1', 'LASSO', 'Ridge', 'OLS', 'Oracle')
colnames(Specificity) = c('beta.vec', 'p', 'N', 'SCAD1', 'LASSO', 'Ridge', 'OLS', 'Oracle')
colnames(Accuracry) = c('beta.vec', 'p', 'N', 'SCAD1', 'LASSO', 'Ridge', 'OLS', 'Oracle')
colnames(MSE) = c('beta.vec', 'p', 'N', 'SCAD1', 'LASSO', 'Ridge', 'OLS', 'Oracle')
#beta.vec : 1 = (1/1, 1/2, 1/3, 1/4, 1/5, 1/6, ... , 1/p) (length : p)
#beta.vec : 2 = (1/1, 1/2, 1/3, 1/4, 1/5, 0, 0, ... ,0) (length : p)
Result Table
## $Sensitivity
## beta.vec p N SCAD1 LASSO Ridge OLS Oracle
## [1,] 1 10 50 0.7150 0.6170 1 1 1
## [2,] 1 10 100 0.8550 0.7360 1 1 1
## [3,] 1 10 200 0.9310 0.8420 1 1 1
## [4,] 1 10 400 0.9750 0.9050 1 1 1
## [5,] 1 10 800 0.9910 0.9640 1 1 1
## [6,] 1 10 1600 0.9980 0.9900 1 1 1
## [7,] 2 10 50 0.7900 0.7600 1 1 1
## [8,] 2 10 100 0.8960 0.8720 1 1 1
## [9,] 2 10 200 0.9620 0.9400 1 1 1
## [10,] 2 10 400 0.9880 0.9880 1 1 1
## [11,] 2 10 800 1.0000 0.9940 1 1 1
## [12,] 2 10 1600 1.0000 1.0000 1 1 1
## [13,] 1 100 50 0.1389 0.0915 1 NA NA
## [14,] 1 100 100 0.2228 0.1299 1 NA NA
## [15,] 1 100 200 0.3111 0.1873 1 1 1
## [16,] 1 100 400 0.4245 0.2682 1 1 1
## [17,] 1 100 800 0.5407 0.3701 1 1 1
## [18,] 1 100 1600 0.6409 0.4736 1 1 1
## [19,] 2 100 50 0.6920 0.6740 1 1 1
## [20,] 2 100 100 0.8280 0.8120 1 1 1
## [21,] 2 100 200 0.8960 0.9160 1 1 1
## [22,] 2 100 400 0.9480 0.9640 1 1 1
## [23,] 2 100 800 0.9940 0.9960 1 1 1
## [24,] 2 100 1600 1.0000 1.0000 1 1 1
##
## $Specificity
## beta.vec p N SCAD1 LASSO Ridge OLS Oracle
## [1,] 1 10 50 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [2,] 1 10 100 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [3,] 1 10 200 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [4,] 1 10 400 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [5,] 1 10 800 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [6,] 1 10 1600 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [7,] 2 10 50 0.6700000 0.6700000 0.6700000 0.6700000 0.6700000
## [8,] 2 10 100 0.6980000 0.6980000 0.6980000 0.6980000 0.6980000
## [9,] 2 10 200 0.6900000 0.6900000 0.6900000 0.6900000 0.6900000
## [10,] 2 10 400 0.7400000 0.7400000 0.7400000 0.7400000 0.7400000
## [11,] 2 10 800 0.7240000 0.7240000 0.7240000 0.7240000 0.7240000
## [12,] 2 10 1600 0.7360000 0.7360000 0.7360000 0.7360000 0.7360000
## [13,] 1 100 50 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [14,] 1 100 100 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [15,] 1 100 200 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [16,] 1 100 400 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [17,] 1 100 800 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [18,] 1 100 1600 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [19,] 2 100 50 0.9327368 0.9327368 0.9327368 0.9327368 0.9327368
## [20,] 2 100 100 0.9298947 0.9298947 0.9298947 0.9298947 0.9298947
## [21,] 2 100 200 0.9336842 0.9336842 0.9336842 0.9336842 0.9336842
## [22,] 2 100 400 0.9310526 0.9310526 0.9310526 0.9310526 0.9310526
## [23,] 2 100 800 0.9258947 0.9258947 0.9258947 0.9258947 0.9258947
## [24,] 2 100 1600 0.9667368 0.9667368 0.9667368 0.9667368 0.9667368
##
## $Accuracry
## beta.vec p N SCAD1 LASSO Ridge OLS Oracle
## [1,] 1 10 50 0.00 0.00 0.00 0.00 0.00
## [2,] 1 10 100 0.00 0.00 0.00 0.00 0.00
## [3,] 1 10 200 0.00 0.00 0.00 0.00 0.00
## [4,] 1 10 400 0.00 0.00 0.00 0.00 0.00
## [5,] 1 10 800 0.00 0.00 0.00 0.00 0.00
## [6,] 1 10 1600 0.00 0.00 0.00 0.00 0.00
## [7,] 2 10 50 0.03 0.05 0.35 0.35 0.35
## [8,] 2 10 100 0.11 0.12 0.34 0.34 0.34
## [9,] 2 10 200 0.20 0.20 0.32 0.32 0.32
## [10,] 2 10 400 0.35 0.37 0.39 0.39 0.39
## [11,] 2 10 800 0.41 0.41 0.41 0.41 0.41
## [12,] 2 10 1600 0.40 0.40 0.40 0.40 0.40
## [13,] 1 100 50 0.00 0.00 0.00 NA NA
## [14,] 1 100 100 0.00 0.00 0.00 NA NA
## [15,] 1 100 200 0.00 0.00 0.00 0.00 0.00
## [16,] 1 100 400 0.00 0.00 0.00 0.00 0.00
## [17,] 1 100 800 0.00 0.00 0.00 0.00 0.00
## [18,] 1 100 1600 0.00 0.00 0.00 0.00 0.00
## [19,] 2 100 50 0.00 0.01 0.14 0.14 0.14
## [20,] 2 100 100 0.03 0.01 0.18 0.18 0.18
## [21,] 2 100 200 0.03 0.08 0.14 0.14 0.14
## [22,] 2 100 400 0.07 0.18 0.23 0.23 0.23
## [23,] 2 100 800 0.19 0.20 0.22 0.22 0.22
## [24,] 2 100 1600 0.50 0.50 0.50 0.50 0.50
##
## $MSE
## beta.vec p N SCAD1 LASSO Ridge OLS
## [1,] 1 10 50 0.404755721 0.29762149 0.30430035 0.424769261
## [2,] 1 10 100 0.193546943 0.16916916 0.18325033 0.175120056
## [3,] 1 10 200 0.089952043 0.10909246 0.12215801 0.082396667
## [4,] 1 10 400 0.043985186 0.06902261 0.07937598 0.041265714
## [5,] 1 10 800 0.020183887 0.03835355 0.05348549 0.019514787
## [6,] 1 10 1600 0.010516975 0.02416221 0.03543160 0.010300534
## [7,] 2 10 50 0.379695780 0.29412156 0.33923435 0.455833338
## [8,] 2 10 100 0.152648753 0.15268398 0.19370512 0.173806923
## [9,] 2 10 200 0.074851364 0.08835470 0.12101552 0.083108538
## [10,] 2 10 400 0.032602428 0.05359687 0.08119414 0.039260311
## [11,] 2 10 800 0.015783755 0.03522195 0.05584505 0.021082927
## [12,] 2 10 1600 0.007187521 0.02189435 0.03829239 0.009915126
## [13,] 1 100 50 0.858364352 0.57109642 1.46219261 NA
## [14,] 1 100 100 0.506657947 0.33572371 1.27083168 NA
## [15,] 1 100 200 0.292157362 0.20298776 0.41710227 1.732178813
## [16,] 1 100 400 0.174489669 0.12942467 0.25153446 0.560721558
## [17,] 1 100 800 0.105209015 0.08314111 0.15621786 0.237567588
## [18,] 1 100 1600 0.066095164 0.05451755 0.09254334 0.109479994
## [19,] 2 100 50 0.454285184 0.36384573 1.35332899 NA
## [20,] 2 100 100 0.244324667 0.19326226 1.21465425 NA
## [21,] 2 100 200 0.119190962 0.11056532 0.42073120 1.656540995
## [22,] 2 100 400 0.061414349 0.06591561 0.25604243 0.551169599
## [23,] 2 100 800 0.022950453 0.03743628 0.15519500 0.243128132
## [24,] 2 100 1600 0.008058934 0.02580570 0.09667600 0.111283250
## Oracle
## [1,] 0.424769261
## [2,] 0.175120056
## [3,] 0.082396667
## [4,] 0.041265714
## [5,] 0.019514787
## [6,] 0.010300534
## [7,] 0.185256653
## [8,] 0.084071776
## [9,] 0.041600756
## [10,] 0.018962789
## [11,] 0.009840577
## [12,] 0.005033200
## [13,] NA
## [14,] NA
## [15,] 1.732178813
## [16,] 0.560721558
## [17,] 0.237567588
## [18,] 0.109479994
## [19,] 0.193336902
## [20,] 0.088951631
## [21,] 0.040252094
## [22,] 0.019302247
## [23,] 0.009619731
## [24,] 0.005470244