## Kelsi Perttula
## HW #2, Problem #3 install.packages('mvtnorm') install.packages('glmnet')
## install.packages('locfit') install.packages('xtable')
## install.packages('bigmemory') install.packages('biganalytics')
## install.packages('plyr') install.packages('multicore')
## install.packages('snowfall') install.packages('rlecuyer')
## install.packages('latticeExtra') install.packages('lasso2')
library(mvtnorm)
library(glmnet)
## Loading required package: Matrix
## Warning: package 'Matrix' was built under R version 3.0.2
## Loaded glmnet 1.9-5
library(locfit)
## locfit 1.5-9.1 2013-03-22
library(xtable)
library(bigmemory)
## Warning: package 'bigmemory' was built under R version 3.0.2
## Loading required package: bigmemory.sri
## Loading required package: BH
##
## bigmemory >= 4.0 is a major revision since 3.1.2; please see packages
## biganalytics and and bigtabulate and http://www.bigmemory.org for more information.
library(biganalytics)
library(plyr)
library(parallel)
library(multicore)
##
## Attaching package: 'multicore'
##
## The following object is masked from 'package:parallel':
##
## mclapply, mcparallel, pvec
library(snowfall)
## Loading required package: snow
##
## Attaching package: 'snow'
##
## The following objects are masked from 'package:parallel':
##
## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
## clusterExport, clusterMap, clusterSplit, makeCluster,
## parApply, parCapply, parLapply, parRapply, parSapply,
## splitIndices, stopCluster
library(rlecuyer)
library(latticeExtra)
## Loading required package: RColorBrewer
## Loading required package: lattice
##
## Attaching package: 'lattice'
##
## The following object is masked from 'package:multicore':
##
## parallel
library(lasso2)
## Warning: package 'lasso2' was built under R version 3.0.2
## R Package to solve regression problems while imposing
## an L1 constraint on the parameters. Based on S-plus Release 2.1
## Copyright (C) 1998, 1999
## Justin Lokhorst <jlokhors@stats.adelaide.edu.au>
## Berwin A. Turlach <bturlach@stats.adelaide.edu.au>
## Bill Venables <wvenable@stats.adelaide.edu.au>
##
## Copyright (C) 2002
## Martin Maechler <maechler@stat.math.ethz.ch>
##
## Attaching package: 'lasso2'
##
## The following object is masked from 'package:locfit':
##
## gcv
library(gplots)
## Loading required package: gtools
## Warning: package 'gtools' was built under R version 3.0.2
## Loading required package: gdata
## gdata: read.xls support for 'XLS' (Excel 97-2004) files ENABLED.
##
## gdata: read.xls support for 'XLSX' (Excel 2007+) files ENABLED.
##
## Attaching package: 'gdata'
##
## The following object is masked from 'package:stats':
##
## nobs
##
## The following object is masked from 'package:utils':
##
## object.size
##
## Loading required package: caTools
## Warning: package 'caTools' was built under R version 3.0.2
## Loading required package: grid
## Loading required package: KernSmooth
## KernSmooth 2.23 loaded
## Copyright M. P. Wand 1997-2009
## Loading required package: MASS
##
## Attaching package: 'gplots'
##
## The following object is masked from 'package:stats':
##
## lowess
# 3a
J <- 10
betas <- c(rep(-J/2, 5) + c(1:5), rep(J/2, 5) - c(5:1))/10
sigma <- 2
rho <- 0.5
nLearn <- 100
nTest <- 1000
set.seed(700)
mu <- rep(0, J)
mu
## [1] 0 0 0 0 0 0 0 0 0 0
sigma <- matrix(1, J, J)
head(sigma)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 1 1 1 1 1 1 1 1 1 1
## [2,] 1 1 1 1 1 1 1 1 1 1
## [3,] 1 1 1 1 1 1 1 1 1 1
## [4,] 1 1 1 1 1 1 1 1 1 1
## [5,] 1 1 1 1 1 1 1 1 1 1
## [6,] 1 1 1 1 1 1 1 1 1 1
for (i in (1:J)) {
for (j in (1:J)) {
sigma[i, j] <- rho^abs(i - j)
}
} #covariance matrix- symmetric
head(sigma)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,] 1.00000 0.5000 0.250 0.125 0.0625 0.03125 0.01562 0.007812 0.003906
## [2,] 0.50000 1.0000 0.500 0.250 0.1250 0.06250 0.03125 0.015625 0.007812
## [3,] 0.25000 0.5000 1.000 0.500 0.2500 0.12500 0.06250 0.031250 0.015625
## [4,] 0.12500 0.2500 0.500 1.000 0.5000 0.25000 0.12500 0.062500 0.031250
## [5,] 0.06250 0.1250 0.250 0.500 1.0000 0.50000 0.25000 0.125000 0.062500
## [6,] 0.03125 0.0625 0.125 0.250 0.5000 1.00000 0.50000 0.250000 0.125000
## [,10]
## [1,] 0.001953
## [2,] 0.003906
## [3,] 0.007812
## [4,] 0.015625
## [5,] 0.031250
## [6,] 0.062500
dim(sigma)
## [1] 10 10
n <- (nLearn + nTest)
X <- rmvnorm(n, mean = mu, sigma = sigma)
X[1:5, 1:5]
## [,1] [,2] [,3] [,4] [,5]
## [1,] -0.1994 -0.3738 -0.4626 -1.4052 0.53602
## [2,] 0.3557 0.1871 0.3822 -1.4330 -1.55267
## [3,] 0.9642 0.7964 1.6387 0.5436 0.62349
## [4,] -0.9820 -0.2320 -1.2188 0.1012 -0.20201
## [5,] 1.8482 -0.9791 0.8889 0.3203 -0.04461
dim(X)
## [1] 1100 10
# Predict Y
Xi_x_Betai <- apply(X, 1, function(x_i) {
sum(x_i * betas)
})
head(Xi_x_Betai)
## [1] -0.36972 0.02421 -1.51876 2.07140 -0.40277 -0.65760
length(Xi_x_Betai)
## [1] 1100
Y <- rnorm(n, mean = Xi_x_Betai, sd = sigma)
head(Y)
## [1] 0.3602 0.3566 -1.7252 1.9168 -0.4603 -0.6117
length(Y)
## [1] 1100
trainSet <- c(rep(1, 100), rep(0, 1000))
simulData <- data.frame(cbind(X, Y, trainSet))
head(simulData)
## V1 V2 V3 V4 V5 V6 V7 V8 V9
## 1 -0.1994 -0.3738 -0.4626 -1.4052 0.53602 0.1112 1.0532 -0.7550 -1.04458
## 2 0.3557 0.1871 0.3822 -1.4330 -1.55267 -0.9714 0.3512 -0.4772 -0.13564
## 3 0.9642 0.7964 1.6387 0.5436 0.62349 0.5082 0.2413 0.6524 -0.09502
## 4 -0.9820 -0.2320 -1.2188 0.1012 -0.20201 0.2792 1.9281 0.9450 0.95381
## 5 1.8482 -0.9791 0.8889 0.3203 -0.04461 1.0019 1.8212 0.1477 -0.03880
## 6 1.4840 -0.7143 0.3336 0.5575 1.14956 0.6682 0.1380 0.5520 0.23514
## V10 Y trainSet
## 1 -1.0890 0.3602 1
## 2 0.6419 0.3566 1
## 3 -1.5954 -1.7252 1
## 4 1.7685 1.9168 1
## 5 0.1314 -0.4603 1
## 6 -0.8764 -0.6117 1
tail(simulData)
## V1 V2 V3 V4 V5 V6 V7 V8
## 1095 1.4900 0.5200 0.03267 1.0320 1.0569 0.4424 -0.7328 0.95312
## 1096 0.8396 0.9081 -0.43137 -1.0270 0.6206 1.2421 0.6608 -0.11600
## 1097 1.4568 2.9236 1.22431 2.0260 -0.5745 -0.9671 0.2199 0.03998
## 1098 -0.7022 -0.5907 1.19678 1.2652 -0.3496 -0.7901 -1.3297 0.07258
## 1099 0.4187 0.1885 -0.95416 -1.0787 0.2849 -0.6431 -0.6361 -1.41249
## 1100 0.3093 1.9245 0.31242 -0.2138 -0.7300 -0.8698 0.9203 0.53280
## V9 V10 Y trainSet
## 1095 -1.1040 0.2557 -0.9696666 0
## 1096 -0.2476 0.7187 -0.1749723 0
## 1097 0.1579 1.2600 -1.5008004 0
## 1098 1.3467 0.7213 0.4471258 0
## 1099 -1.7466 -0.5759 -1.5492255 0
## 1100 -0.1381 0.3905 -0.0002092 0
dim(simulData)
## [1] 1100 12
colnames(simulData) <- c(paste("X", c(1:10)), "Y", "train")
learnSet <- simulData[c(1:100), ]
head(learnSet)
## X 1 X 2 X 3 X 4 X 5 X 6 X 7 X 8 X 9
## 1 -0.1994 -0.3738 -0.4626 -1.4052 0.53602 0.1112 1.0532 -0.7550 -1.04458
## 2 0.3557 0.1871 0.3822 -1.4330 -1.55267 -0.9714 0.3512 -0.4772 -0.13564
## 3 0.9642 0.7964 1.6387 0.5436 0.62349 0.5082 0.2413 0.6524 -0.09502
## 4 -0.9820 -0.2320 -1.2188 0.1012 -0.20201 0.2792 1.9281 0.9450 0.95381
## 5 1.8482 -0.9791 0.8889 0.3203 -0.04461 1.0019 1.8212 0.1477 -0.03880
## 6 1.4840 -0.7143 0.3336 0.5575 1.14956 0.6682 0.1380 0.5520 0.23514
## X 10 Y train
## 1 -1.0890 0.3602 1
## 2 0.6419 0.3566 1
## 3 -1.5954 -1.7252 1
## 4 1.7685 1.9168 1
## 5 0.1314 -0.4603 1
## 6 -0.8764 -0.6117 1
dim(learnSet)
## [1] 100 12
testSet <- simulData[c(101:1100), ]
head(testSet)
## X 1 X 2 X 3 X 4 X 5 X 6 X 7 X 8
## 101 -0.2169 -1.03690 -1.4390 -0.09723 -0.82080 -1.2752 -0.6306 -1.12867
## 102 -0.6578 -0.08393 0.6663 0.77206 0.35673 1.2715 0.5064 0.01925
## 103 -0.6823 0.38979 0.1094 -0.28733 0.87111 1.2303 1.7231 0.59997
## 104 -0.0114 0.91976 2.1094 0.45209 0.36691 -0.5161 -0.9519 -1.45527
## 105 1.2754 -1.38717 -0.9441 -1.34031 -1.02475 0.0624 -1.9434 -1.37122
## 106 2.0168 1.06797 0.3548 -0.08322 -0.04206 0.7172 -0.4549 0.52513
## X 9 X 10 Y train
## 101 -1.71410 -2.1959 -1.7988 0
## 102 -0.08939 0.8454 0.3121 0
## 103 1.12140 -0.3424 0.8850 0
## 104 -1.93791 -2.7194 -2.9143 0
## 105 -2.24218 -0.7091 -1.1787 0
## 106 0.61785 0.6410 -0.6834 0
dim(testSet)
## [1] 1000 12
# EDA of learningSet
pairs(learnSet[, 1:11], col = c(5), main = "Simulated Data-Pairs - Learning Set")
myPalette <- colorRampPalette(c("blue", "green"), space = "rgb")
heatmap.2(cor(learnSet[, 1:11]), col = myPalette, main = "Learning Set")
# Outside of diagonal, correlation increases away from diagonal
par(mfrow = c(2, 5))
for (j in 1:10) {
plot(learnSet[, j], learnSet[, 11], col = c(7, 3)[learnSet$train * 1 + 1],
xlab = names(learnSet)[j], ylab = names(learnSet)[11], main = "", sub = paste("Corr=",
round(cor(learnSet[, c(j, 11)])[1, 2], 2), sep = ""), cex.main = 1.1,
cex.lab = 0.9)
lines(lowess(learnSet[, 11] ~ learnSet[, j]), col = 1)
abline(lm(learnSet[, 11] ~ learnSet[, j])$coef, col = 6)
}
par(mfrow = c(1, 1))
# Elastic Net
lambdas <- c(0:100)
revLambdas <- rev(lambdas)
learnCovars <- as.matrix(learnSet[, c(1:10)])
dim(learnCovars)
## [1] 100 10
learnY <- learnSet$Y
# a- Ridge
alphaRidge <- 0
lambdasRidge <- lambdas/(100)
lambdasRidge
## [1] 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 0.11 0.12 0.13
## [15] 0.14 0.15 0.16 0.17 0.18 0.19 0.20 0.21 0.22 0.23 0.24 0.25 0.26 0.27
## [29] 0.28 0.29 0.30 0.31 0.32 0.33 0.34 0.35 0.36 0.37 0.38 0.39 0.40 0.41
## [43] 0.42 0.43 0.44 0.45 0.46 0.47 0.48 0.49 0.50 0.51 0.52 0.53 0.54 0.55
## [57] 0.56 0.57 0.58 0.59 0.60 0.61 0.62 0.63 0.64 0.65 0.66 0.67 0.68 0.69
## [71] 0.70 0.71 0.72 0.73 0.74 0.75 0.76 0.77 0.78 0.79 0.80 0.81 0.82 0.83
## [85] 0.84 0.85 0.86 0.87 0.88 0.89 0.90 0.91 0.92 0.93 0.94 0.95 0.96 0.97
## [99] 0.98 0.99 1.00
fitRidge <- glmnet(learnCovars, learnY, standardize = F, family = "gaussian",
intercept = F, lambda = lambdasRidge, alpha = alphaRidge)
(fitRidge$beta[, c(length(lambdas):1)])[1:5, 1:5]
## 5 x 5 sparse Matrix of class "dgCMatrix"
## s100 s99 s98 s97 s96
## X 1 -0.370552 -0.366922 -0.363391 -0.35995 -0.35661
## X 2 -0.331785 -0.331307 -0.330729 -0.33006 -0.32932
## X 3 -0.172706 -0.171321 -0.170003 -0.16875 -0.16754
## X 4 -0.100484 -0.100043 -0.099651 -0.09930 -0.09899
## X 5 -0.006436 -0.008207 -0.009844 -0.01136 -0.01278
alphaLasso <- 1
lambdasLasso <- lambdas/(2 * 100)
fitLasso <- glmnet(x = learnCovars, y = learnY, standardize = F, family = "gaussian",
lambda = lambdasLasso, alpha = alphaLasso)
fitLasso$beta[, c(length(lambdas):1)]
## 10 x 101 sparse Matrix of class "dgCMatrix"
## [[ suppressing 101 column names 's100', 's99', 's98' ... ]]
##
## X 1 -0.373937 -0.36641 -0.35937 -0.35234 -0.34530 -0.33826 -0.331222
## X 2 -0.326159 -0.32797 -0.32888 -0.32979 -0.33070 -0.33161 -0.332527
## X 3 -0.175653 -0.17184 -0.16848 -0.16511 -0.16174 -0.15838 -0.155008
## X 4 -0.102547 -0.09491 -0.09055 -0.08619 -0.08184 -0.07748 -0.073120
## X 5 0.007185 . . . . . .
## X 6 0.046104 0.04334 0.03637 0.02941 0.02244 0.01548 0.008512
## X 7 0.093623 0.09467 0.09535 0.09602 0.09670 0.09738 0.098051
## X 8 0.220816 0.21714 0.21512 0.21311 0.21109 0.20908 0.207062
## X 9 0.208335 0.20493 0.20083 0.19672 0.19261 0.18851 0.184403
## X 10 0.448573 0.44452 0.44170 0.43888 0.43606 0.43324 0.430415
##
## X 1 -0.324167 -0.31890 -0.31407 -0.30925 -0.30442 -0.29960 -0.29478
## X 2 -0.333528 -0.33394 -0.33420 -0.33446 -0.33472 -0.33498 -0.33523
## X 3 -0.151598 -0.14808 -0.14459 -0.14111 -0.13762 -0.13414 -0.13066
## X 4 -0.068779 -0.06579 -0.06315 -0.06051 -0.05788 -0.05524 -0.05260
## X 5 . . . . . . .
## X 6 0.001574 . . . . . .
## X 7 0.098662 0.09710 0.09498 0.09287 0.09075 0.08863 0.08651
## X 8 0.205032 0.20283 0.20054 0.19826 0.19597 0.19369 0.19140
## X 9 0.180448 0.17684 0.17328 0.16972 0.16615 0.16259 0.15903
## X 10 0.427511 0.42361 0.41946 0.41530 0.41115 0.40699 0.40283
##
## X 1 -0.28995 -0.28513 -0.28030 -0.27548 -0.27065 -0.26583 -0.26101
## X 2 -0.33549 -0.33575 -0.33601 -0.33627 -0.33653 -0.33679 -0.33704
## X 3 -0.12717 -0.12369 -0.12020 -0.11672 -0.11323 -0.10975 -0.10626
## X 4 -0.04996 -0.04733 -0.04469 -0.04205 -0.03941 -0.03678 -0.03414
## X 5 . . . . . . .
## X 6 . . . . . . .
## X 7 0.08439 0.08227 0.08015 0.07803 0.07591 0.07380 0.07168
## X 8 0.18912 0.18683 0.18455 0.18226 0.17998 0.17769 0.17541
## X 9 0.15547 0.15190 0.14834 0.14478 0.14122 0.13765 0.13409
## X 10 0.39868 0.39452 0.39036 0.38621 0.38205 0.37790 0.37374
##
## X 1 -0.25618 -0.25136 -0.24653 -0.24171 -0.23689 -0.23206 -0.22724
## X 2 -0.33730 -0.33756 -0.33782 -0.33808 -0.33834 -0.33860 -0.33885
## X 3 -0.10278 -0.09929 -0.09581 -0.09232 -0.08884 -0.08536 -0.08187
## X 4 -0.03150 -0.02886 -0.02623 -0.02359 -0.02095 -0.01832 -0.01568
## X 5 . . . . . . .
## X 6 . . . . . . .
## X 7 0.06956 0.06744 0.06532 0.06320 0.06108 0.05896 0.05685
## X 8 0.17312 0.17084 0.16856 0.16627 0.16399 0.16170 0.15942
## X 9 0.13053 0.12697 0.12340 0.11984 0.11628 0.11272 0.10916
## X 10 0.36958 0.36543 0.36127 0.35712 0.35296 0.34880 0.34465
##
## X 1 -0.22241 -0.21759 -0.212764 -0.207943 -0.203122 -0.19826 -0.19326
## X 2 -0.33911 -0.33937 -0.339630 -0.339886 -0.340057 -0.34044 -0.34021
## X 3 -0.07839 -0.07490 -0.071417 -0.067934 -0.064526 -0.06088 -0.05637
## X 4 -0.01304 -0.01040 -0.007765 -0.005128 -0.002465 . .
## X 5 . . . . . . .
## X 6 . . . . . . .
## X 7 0.05473 0.05261 0.050489 0.048369 0.046355 0.04405 0.04175
## X 8 0.15713 0.15485 0.152561 0.150278 0.147955 0.14574 0.14392
## X 9 0.10559 0.10203 0.098468 0.094903 0.091220 0.08789 0.08445
## X 10 0.34049 0.33633 0.332178 0.328023 0.323939 0.31965 0.31526
##
## X 1 -0.18827 -0.18327 -0.17828 -0.17328 -0.16829 -0.16329 -0.15830
## X 2 -0.33998 -0.33975 -0.33953 -0.33930 -0.33907 -0.33884 -0.33861
## X 3 -0.05185 -0.04734 -0.04283 -0.03831 -0.03380 -0.02929 -0.02478
## X 4 . . . . . . .
## X 5 . . . . . . .
## X 6 . . . . . . .
## X 7 0.03945 0.03716 0.03486 0.03256 0.03026 0.02797 0.02567
## X 8 0.14210 0.14029 0.13847 0.13665 0.13483 0.13302 0.13120
## X 9 0.08101 0.07758 0.07414 0.07070 0.06726 0.06383 0.06039
## X 10 0.31087 0.30649 0.30210 0.29771 0.29333 0.28894 0.28455
##
## X 1 -0.15330 -0.14830 -0.14331 -0.138316 -0.133279 -0.12861 -0.12407
## X 2 -0.33839 -0.33816 -0.33793 -0.337698 -0.337466 -0.33589 -0.33311
## X 3 -0.02026 -0.01575 -0.01124 -0.006728 -0.002231 . .
## X 4 . . . . . . .
## X 5 . . . . . . .
## X 6 . . . . . . .
## X 7 0.02337 0.02108 0.01878 0.016480 0.014268 0.01225 0.01063
## X 8 0.12938 0.12756 0.12575 0.123930 0.122091 0.12013 0.11793
## X 9 0.05695 0.05351 0.05007 0.046632 0.043104 0.03946 0.03557
## X 10 0.28017 0.27578 0.27139 0.267009 0.262663 0.25874 0.25529
##
## X 1 -0.119537 -0.115001 -0.110464 -0.10593 -0.101391 -0.0968219 -0.092283
## X 2 -0.330334 -0.327555 -0.324777 -0.32200 -0.319221 -0.3164633 -0.313607
## X 3 . . . . . . .
## X 4 . . . . . . .
## X 5 . . . . . . .
## X 6 . . . . . . .
## X 7 0.009021 0.007407 0.005794 0.00418 0.002566 0.0009113 .
## X 8 0.115741 0.113548 0.111356 0.10916 0.106970 0.1047601 0.102357
## X 9 0.031689 0.027806 0.023922 0.02004 0.016156 0.0122919 0.008043
## X 10 0.251831 0.248376 0.244920 0.24146 0.238009 0.2345555 0.231215
##
## X 1 -0.087638 -0.08303 -0.07863 -0.07423 -0.06983 -0.06544 -0.06104
## X 2 -0.310784 -0.30810 -0.30588 -0.30367 -0.30145 -0.29924 -0.29703
## X 3 . . . . . . .
## X 4 . . . . . . .
## X 5 . . . . . . .
## X 6 . . . . . . .
## X 7 . . . . . . .
## X 8 0.099510 0.09636 0.09206 0.08777 0.08347 0.07917 0.07488
## X 9 0.003561 . . . . . .
## X 10 0.227934 0.22423 0.21889 0.21355 0.20821 0.20287 0.19753
##
## X 1 -0.05664 -0.05224 -0.04784 -0.04344 -0.03904 -0.03465 -0.03025
## X 2 -0.29481 -0.29260 -0.29038 -0.28817 -0.28596 -0.28374 -0.28153
## X 3 . . . . . . .
## X 4 . . . . . . .
## X 5 . . . . . . .
## X 6 . . . . . . .
## X 7 . . . . . . .
## X 8 0.07058 0.06629 0.06199 0.05769 0.05340 0.04910 0.04480
## X 9 . . . . . . .
## X 10 0.19219 0.18685 0.18151 0.17617 0.17082 0.16548 0.16014
##
## X 1 -0.02585 -0.02145 -0.01705 -0.01265 -0.008255 -0.003858 .
## X 2 -0.27931 -0.27710 -0.27488 -0.27267 -0.270457 -0.268242 -0.26579
## X 3 . . . . . . .
## X 4 . . . . . . .
## X 5 . . . . . . .
## X 6 . . . . . . .
## X 7 . . . . . . .
## X 8 0.04051 0.03621 0.03192 0.02762 0.023323 0.019028 0.01478
## X 9 . . . . . . .
## X 10 0.15480 0.14946 0.14412 0.13878 0.133440 0.128099 0.12274
##
## X 1 . . . . . . .
## X 2 -0.26166 -0.2575 -0.253385 -0.2492 -0.24504 -0.24083 -0.23663
## X 3 . . . . . . .
## X 4 . . . . . . .
## X 5 . . . . . . .
## X 6 . . . . . . .
## X 7 . . . . . . .
## X 8 0.01089 0.0070 0.003101 . . . .
## X 9 . . . . . . .
## X 10 0.11722 0.1117 0.106182 0.1004 0.09388 0.08732 0.08075
##
## X 1 . . . . . . .
## X 2 -0.23242 -0.22822 -0.22402 -0.21981 -0.21561 -0.21140 -0.2072
## X 3 . . . . . . .
## X 4 . . . . . . .
## X 5 . . . . . . .
## X 6 . . . . . . .
## X 7 . . . . . . .
## X 8 . . . . . . .
## X 9 . . . . . . .
## X 10 0.07419 0.06762 0.06106 0.05449 0.04793 0.04136 0.0348
##
## X 1 . . . . . . . .
## X 2 -0.20300 -0.19879 -0.1946 -0.190383 -0.186179 -0.1815 -0.1767 -0.1719
## X 3 . . . . . . . .
## X 4 . . . . . . . .
## X 5 . . . . . . . .
## X 6 . . . . . . . .
## X 7 . . . . . . . .
## X 8 . . . . . . . .
## X 9 . . . . . . . .
## X 10 0.02823 0.02167 0.0151 0.008539 0.001974 . . .
##
## X 1 . .
## X 2 -0.167 -0.1622
## X 3 . .
## X 4 . .
## X 5 . .
## X 6 . .
## X 7 . .
## X 8 . .
## X 9 . .
## X 10 . .
fitLasso$lambda
## [1] 0.500 0.495 0.490 0.485 0.480 0.475 0.470 0.465 0.460 0.455 0.450
## [12] 0.445 0.440 0.435 0.430 0.425 0.420 0.415 0.410 0.405 0.400 0.395
## [23] 0.390 0.385 0.380 0.375 0.370 0.365 0.360 0.355 0.350 0.345 0.340
## [34] 0.335 0.330 0.325 0.320 0.315 0.310 0.305 0.300 0.295 0.290 0.285
## [45] 0.280 0.275 0.270 0.265 0.260 0.255 0.250 0.245 0.240 0.235 0.230
## [56] 0.225 0.220 0.215 0.210 0.205 0.200 0.195 0.190 0.185 0.180 0.175
## [67] 0.170 0.165 0.160 0.155 0.150 0.145 0.140 0.135 0.130 0.125 0.120
## [78] 0.115 0.110 0.105 0.100 0.095 0.090 0.085 0.080 0.075 0.070 0.065
## [89] 0.060 0.055 0.050 0.045 0.040 0.035 0.030 0.025 0.020 0.015 0.010
## [100] 0.005 0.000
# Elastic Net
alphaEnet <- 1/3
lambdasEnet <- 3 * (lambdas)/(2 * 100)
fitEnet <- glmnet(x = learnCovars, y = learnY, standardize = F, family = "gaussian",
lambda = lambdasEnet, alpha = alphaEnet)
fitEnet$beta[, c(length(lambdas):1)]
## 10 x 101 sparse Matrix of class "dgCMatrix"
## [[ suppressing 101 column names 's100', 's99', 's98' ... ]]
##
## X 1 -0.373937 -0.36295 -0.35286 -0.34315 -0.33378 -0.32473 -0.315980
## X 2 -0.326159 -0.32735 -0.32737 -0.32715 -0.32671 -0.32607 -0.325243
## X 3 -0.175653 -0.17049 -0.16592 -0.16148 -0.15716 -0.15295 -0.148857
## X 4 -0.102547 -0.09541 -0.09167 -0.08803 -0.08449 -0.08103 -0.077655
## X 5 0.007185 . . . . . .
## X 6 0.046104 0.04283 0.03569 0.02883 0.02223 0.01589 0.009781
## X 7 0.093623 0.09535 0.09650 0.09746 0.09826 0.09890 0.099396
## X 8 0.220816 0.21579 0.21251 0.20929 0.20613 0.20302 0.199962
## X 9 0.208335 0.20484 0.20069 0.19658 0.19251 0.18847 0.184482
## X 10 0.448573 0.43897 0.43078 0.42277 0.41492 0.40722 0.399689
##
## X 1 -0.307446 -0.29984 -0.29341 -0.28712 -0.28096 -0.27494 -0.26904
## X 2 -0.324293 -0.32301 -0.32133 -0.31956 -0.31771 -0.31579 -0.31380
## X 3 -0.144909 -0.14090 -0.13707 -0.13334 -0.12969 -0.12613 -0.12265
## X 4 -0.074314 -0.07156 -0.06961 -0.06765 -0.06569 -0.06373 -0.06178
## X 5 . . . . . . .
## X 6 0.003842 . . . . . .
## X 7 0.099823 0.09928 0.09747 0.09566 0.09386 0.09206 0.09025
## X 8 0.196926 0.19389 0.19066 0.18745 0.18428 0.18114 0.17802
## X 9 0.180551 0.17681 0.17311 0.16942 0.16573 0.16205 0.15838
## X 10 0.392303 0.38473 0.37680 0.36906 0.36151 0.35414 0.34694
##
## X 1 -0.26327 -0.25761 -0.25207 -0.24663 -0.24130 -0.23606 -0.23092
## X 2 -0.31176 -0.30965 -0.30750 -0.30531 -0.30307 -0.30080 -0.29850
## X 3 -0.11925 -0.11591 -0.11265 -0.10944 -0.10630 -0.10321 -0.10018
## X 4 -0.05982 -0.05786 -0.05591 -0.05396 -0.05202 -0.05008 -0.04814
## X 5 . . . . . . .
## X 6 . . . . . . .
## X 7 0.08845 0.08665 0.08485 0.08304 0.08124 0.07944 0.07764
## X 8 0.17494 0.17189 0.16886 0.16586 0.16289 0.15994 0.15702
## X 9 0.15473 0.15109 0.14747 0.14388 0.14031 0.13676 0.13324
## X 10 0.33990 0.33301 0.32626 0.31966 0.31319 0.30685 0.30064
##
## X 1 -0.22588 -0.22092 -0.21605 -0.21127 -0.20656 -0.20194 -0.19739
## X 2 -0.29616 -0.29381 -0.29143 -0.28902 -0.28661 -0.28417 -0.28172
## X 3 -0.09720 -0.09427 -0.09138 -0.08855 -0.08575 -0.08300 -0.08029
## X 4 -0.04622 -0.04429 -0.04237 -0.04046 -0.03856 -0.03666 -0.03477
## X 5 . . . . . . .
## X 6 . . . . . . .
## X 7 0.07584 0.07404 0.07224 0.07044 0.06864 0.06684 0.06505
## X 8 0.15412 0.15125 0.14840 0.14558 0.14278 0.14000 0.13725
## X 9 0.12974 0.12628 0.12284 0.11943 0.11605 0.11270 0.10938
## X 10 0.29454 0.28856 0.28269 0.27693 0.27127 0.26571 0.26025
##
## X 1 -0.19291 -0.18851 -0.18417 -0.17991 -0.17570 -0.17157 -0.16749
## X 2 -0.27926 -0.27679 -0.27431 -0.27183 -0.26933 -0.26684 -0.26434
## X 3 -0.07762 -0.07498 -0.07239 -0.06982 -0.06729 -0.06479 -0.06233
## X 4 -0.03289 -0.03102 -0.02915 -0.02729 -0.02544 -0.02360 -0.02177
## X 5 . . . . . . .
## X 6 . . . . . . .
## X 7 0.06326 0.06147 0.05968 0.05789 0.05611 0.05433 0.05255
## X 8 0.13451 0.13180 0.12912 0.12645 0.12381 0.12118 0.11858
## X 9 0.10609 0.10283 0.09960 0.09640 0.09323 0.09009 0.08698
## X 10 0.25488 0.24960 0.24441 0.23930 0.23428 0.22933 0.22447
##
## X 1 -0.16347 -0.15951 -0.15561 -0.15178 -0.14799 -0.14425 -0.140559
## X 2 -0.26184 -0.25934 -0.25684 -0.25432 -0.25182 -0.24933 -0.246834
## X 3 -0.05989 -0.05748 -0.05510 -0.05275 -0.05043 -0.04813 -0.045854
## X 4 -0.01995 -0.01813 -0.01632 -0.01453 -0.01274 -0.01096 -0.009193
## X 5 . . . . . . .
## X 6 . . . . . . .
## X 7 0.05078 0.04901 0.04725 0.04546 0.04371 0.04195 0.040208
## X 8 0.11600 0.11343 0.11089 0.10837 0.10587 0.10338 0.100916
## X 9 0.08390 0.08085 0.07782 0.07485 0.07188 0.06895 0.066039
## X 10 0.21968 0.21496 0.21032 0.20574 0.20123 0.19679 0.192415
##
## X 1 -0.13692 -0.133331 -0.129789 -0.126294 -0.1228271 -0.11934 -0.11587
## X 2 -0.24435 -0.241861 -0.239382 -0.236908 -0.2344404 -0.23178 -0.22904
## X 3 -0.04361 -0.041383 -0.039182 -0.037004 -0.0348579 -0.03237 -0.02980
## X 4 -0.00743 -0.005677 -0.003933 -0.002198 -0.0004661 . .
## X 5 . . . . . . .
## X 6 . . . . . . .
## X 7 0.03847 0.036728 0.034996 0.033268 0.0315683 0.02981 0.02808
## X 8 0.09847 0.096041 0.093631 0.091240 0.0888631 0.08663 0.08445
## X 9 0.06316 0.060308 0.057484 0.054687 0.0518968 0.04922 0.04656
## X 10 0.18810 0.183846 0.179652 0.175517 0.1714430 0.16737 0.16334
##
## X 1 -0.11244 -0.10905 -0.10570 -0.10239 -0.09912 -0.09589 -0.09270
## X 2 -0.22630 -0.22358 -0.22087 -0.21817 -0.21548 -0.21281 -0.21014
## X 3 -0.02726 -0.02474 -0.02226 -0.01980 -0.01736 -0.01495 -0.01257
## X 4 . . . . . . .
## X 5 . . . . . . .
## X 6 . . . . . . .
## X 7 0.02635 0.02463 0.02292 0.02122 0.01952 0.01783 0.01615
## X 8 0.08229 0.08013 0.07799 0.07587 0.07376 0.07166 0.06957
## X 9 0.04392 0.04131 0.03872 0.03615 0.03361 0.03110 0.02861
## X 10 0.15936 0.15544 0.15157 0.14775 0.14399 0.14027 0.13659
##
## X 1 -0.08955 -0.086428 -0.08335 -0.080299 -0.077267 -0.074310 -0.07136
## X 2 -0.20749 -0.204855 -0.20223 -0.199617 -0.197023 -0.194041 -0.19078
## X 3 -0.01021 -0.007871 -0.00556 -0.003271 -0.001005 . .
## X 4 . . . . . . .
## X 5 . . . . . . .
## X 6 . . . . . . .
## X 7 0.01447 0.012807 0.01115 0.009496 0.007863 0.006301 0.00482
## X 8 0.06750 0.065441 0.06339 0.061361 0.059339 0.057275 0.05518
## X 9 0.02614 0.023690 0.02127 0.018866 0.016475 0.014135 0.01182
## X 10 0.13297 0.129390 0.12586 0.122364 0.118918 0.115597 0.11239
##
## X 1 -0.068442 -0.065553 -0.0626892 -0.059845 -5.701e-02 -0.05433 -0.05167
## X 2 -0.187543 -0.184337 -0.1811580 -0.177987 -1.748e-01 -0.17189 -0.16897
## X 3 . . . . . . .
## X 4 . . . . . . .
## X 5 . . . . . . .
## X 6 . . . . . . .
## X 7 0.003342 0.001867 0.0003943 . . . .
## X 8 0.053099 0.051034 0.0489832 0.046665 4.426e-02 0.04133 0.03844
## X 9 0.009518 0.007243 0.0049798 0.002523 6.573e-06 . .
## X 10 0.109209 0.106068 0.1029637 0.099879 9.682e-02 0.09319 0.08960
##
## X 1 -0.04904 -0.04643 -0.04385 -0.04129 -0.03876 -0.03625 -0.03377
## X 2 -0.16607 -0.16320 -0.16034 -0.15751 -0.15470 -0.15192 -0.14915
## X 3 . . . . . . .
## X 4 . . . . . . .
## X 5 . . . . . . .
## X 6 . . . . . . .
## X 7 . . . . . . .
## X 8 0.03557 0.03274 0.02993 0.02715 0.02440 0.02168 0.01898
## X 9 . . . . . . .
## X 10 0.08605 0.08255 0.07909 0.07567 0.07229 0.06896 0.06566
##
## X 1 -0.03130 -0.02886 -0.02645 -0.024053 -0.021680 -0.019329 -0.016998
## X 2 -0.14641 -0.14369 -0.14099 -0.138307 -0.135648 -0.133009 -0.130391
## X 3 . . . . . . .
## X 4 . . . . . . .
## X 5 . . . . . . .
## X 6 . . . . . . .
## X 7 . . . . . . .
## X 8 0.01631 0.01367 0.01105 0.008456 0.005889 0.003346 0.000828
## X 9 . . . . . . .
## X 10 0.06240 0.05918 0.05599 0.052843 0.049731 0.046654 0.043612
##
## X 1 -0.01477 -0.01260 -0.01044 -0.008303 -0.006182 -0.004078 -0.001991
## X 2 -0.12774 -0.12509 -0.12247 -0.119859 -0.117272 -0.114706 -0.112161
## X 3 . . . . . . .
## X 4 . . . . . . .
## X 5 . . . . . . .
## X 6 . . . . . . .
## X 7 . . . . . . .
## X 8 . . . . . . .
## X 9 . . . . . . .
## X 10 0.04040 0.03712 0.03388 0.030679 0.027517 0.024393 0.021306
##
## X 1 . . . . . . .
## X 2 -0.10962 -0.1066 -0.10368 -0.100753 -0.097853 -0.094981 -0.092135
## X 3 . . . . . . .
## X 4 . . . . . . .
## X 5 . . . . . . .
## X 6 . . . . . . .
## X 7 . . . . . . .
## X 8 . . . . . . .
## X 9 . . . . . . .
## X 10 0.01825 0.0152 0.01218 0.009192 0.006241 0.003323 0.000439
##
## X 1 . . .
## X 2 -0.0892 -0.08626 -0.08336
## X 3 . . .
## X 4 . . .
## X 5 . . .
## X 6 . . .
## X 7 . . .
## X 8 . . .
## X 9 . . .
## X 10 . . .
fitEnet$lambda
## [1] 1.500 1.485 1.470 1.455 1.440 1.425 1.410 1.395 1.380 1.365 1.350
## [12] 1.335 1.320 1.305 1.290 1.275 1.260 1.245 1.230 1.215 1.200 1.185
## [23] 1.170 1.155 1.140 1.125 1.110 1.095 1.080 1.065 1.050 1.035 1.020
## [34] 1.005 0.990 0.975 0.960 0.945 0.930 0.915 0.900 0.885 0.870 0.855
## [45] 0.840 0.825 0.810 0.795 0.780 0.765 0.750 0.735 0.720 0.705 0.690
## [56] 0.675 0.660 0.645 0.630 0.615 0.600 0.585 0.570 0.555 0.540 0.525
## [67] 0.510 0.495 0.480 0.465 0.450 0.435 0.420 0.405 0.390 0.375 0.360
## [78] 0.345 0.330 0.315 0.300 0.285 0.270 0.255 0.240 0.225 0.210 0.195
## [89] 0.180 0.165 0.150 0.135 0.120 0.105 0.090 0.075 0.060 0.045 0.030
## [100] 0.015 0.000
# degrees of freedom vs lamdas
ridgeDf <- sapply(lambdas, function(lambda) {
sum(diag(learnCovars %*% solve(t(learnCovars) %*% learnCovars + lambda *
diag(x = 1, J, J)) %*% t(learnCovars)))
})
matplot(lambdas, cbind(ridgeDf, fitLasso$df[c(length(lambdas):1)], fitEnet$df[c(length(lambdas):1)]),
type = "l", lty = 1, col = 4:7, xlab = expression(lambda), ylab = "Degrees of Freedom",
main = "", cex.lab = 0.6, cex.axis = 0.6)
legend("bottomleft", c("Ridge", "Lasso", "Elastic Net"), lty = 1, col = 4:7,
cex = 0.5, box.col = NA)
# DF for elastic net appears between ridge and Lasso as it is a compromise
# of both methods
# 3b
matplot(revLambdas, t(fitRidge$beta), type = "l", lty = 1, col = 1:12, xlab = expression(lambda),
ylab = "Fitted Coefficients - Ridge", main = "", cex.lab = 0.8, cex.axis = 0.8)
legend("top", c(paste0("X", c(1:10))), lty = 1, col = 1:12, cex = 0.6, ncol = 5)
matplot(revLambdas, t(fitLasso$beta), type = "l", lty = 1, col = 1:12, xlab = expression(lambda),
ylab = "Fitted Coefficients - Lasso", main = "", cex.lab = 0.8, cex.axis = 0.8)
legend("top", c(paste0("X", c(1:10))), lty = 1, col = 1:12, cex = 0.6, ncol = 5)
matplot(revLambdas, t(fitEnet$beta), type = "l", lty = 1, col = 1:12, xlab = expression(lambda),
ylab = "Fitted Coefficients - Enet", main = "", cex.lab = 0.8, cex.axis = 0.8)
legend("top", c(paste0("X", c(1:10))), lty = 1, col = 1:12, cex = 0.6, ncol = 5)
# MSE
yHatRidge <- predict(fitRidge, learnCovars, type = "response")
yHatRidge[1:5, 1:5]
## s0 s1 s2 s3 s4
## 1 -0.098843 -0.099877 -0.100871 -0.101879 -0.102902
## 2 0.076973 0.076919 0.076827 0.076730 0.076630
## 3 -0.862281 -0.865536 -0.868894 -0.872282 -0.875700
## 4 1.254352 1.258792 1.263349 1.267940 1.272567
## 5 -0.006145 -0.006608 -0.007129 -0.007662 -0.008206
dim(yHatRidge)
## [1] 100 101
yHatLasso <- predict(fitLasso, learnCovars, type = "response")
yHatLasso[1:5, 1:5]
## s0 s1 s2 s3 s4
## 1 0.092224 9.366e-02 0.095094 0.09653 0.097964
## 2 0.001251 -2.626e-05 -0.001303 -0.00258 -0.003857
## 3 -0.097565 -1.018e-01 -0.106010 -0.11023 -0.114455
## 4 0.069237 6.999e-02 0.070737 0.07149 0.072236
## 5 0.190413 1.948e-01 0.199137 0.20350 0.207861
dim(yHatLasso)
## [1] 100 101
yHatEnet <- predict(fitEnet, learnCovars, type = "response")
yHatEnet[1:5, 1:5]
## s0 s1 s2 s3 s4
## 1 0.06883 0.06969 0.07056 0.07098 0.06885
## 2 0.02207 0.02131 0.02053 0.02006 0.02133
## 3 -0.02871 -0.03125 -0.03381 -0.03705 -0.04397
## 4 0.05701 0.05746 0.05792 0.05918 0.06489
## 5 0.11929 0.12191 0.12456 0.12729 0.13041
dim(yHatEnet)
## [1] 100 101
mseRidge <- (learnY - yHatRidge)^2
mseLasso <- (learnY - yHatLasso)^2
mseEnet <- (learnY - yHatEnet)^2
mseRidgeMean <- apply(mseRidge, 2, mean)
mseLassoMean <- apply(mseLasso, 2, mean)
mseEnetMean <- apply(mseEnet, 2, mean)
matplot(revLambdas, cbind(mseRidgeMean, mseLassoMean, mseEnetMean), type = "l",
lty = 1, col = 4:7, xlab = expression(lambda), ylab = "MSE[Y]", main = "",
cex.lab = 0.8)
legend("topleft", c("Ridge", "Lasso", "Elastic Net"), lty = 1, col = 4:7, cex = 0.7)
# 3c Test Set w MSE
testCovars <- as.matrix(testSet[, c(1:J)])
testY <- testSet$Y
yhatRidgeTest <- predict(fitRidge, testCovars, type = "response")
yhatLassoTest <- predict(fitLasso, testCovars, type = "response")
yhatEnetTest <- predict(fitEnet, testCovars, type = "response")
mseRidgeTest <- (testY - yhatRidgeTest)^2
mseLassoTest <- (testY - yhatLassoTest)^2
mseEnetTest <- (testY - yhatEnetTest)^2
mseRidgeTestMean <- apply(mseRidgeTest, 2, mean)
mseLassoTestMean <- apply(mseLassoTest, 2, mean)
mseEnetTestMean <- apply(mseEnetTest, 2, mean)
lminRidgeTest <- lambdas[which.min(mseRidgeTestMean)]
lminLassoTest <- lambdas[which.min(mseLassoTestMean)]
lminEnetTest <- lambdas[which.min(mseEnetTestMean)]
lminRidgeTest
## [1] 100
lminLassoTest
## [1] 100
lminEnetTest
## [1] 100
matplot(revLambdas, cbind(mseRidgeMean, mseRidgeTestMean, mseLassoMean, mseLassoTestMean,
mseEnetMean, mseEnetTestMean), type = "l", col = rep(2:4, each = 2), lty = rep(2:1,
3), xlab = expression(lambda), ylab = "MSE[Y]", main = "")
legend("topleft", c("Learning Set", "Test Set"), lty = 2:1, cex = 0.5, box.col = NA)
legend("top", c("Ridge", "Lasso", "Elastic Net"), lty = 1, col = 2:4, cex = 0.5,
box.col = NA)
# 3d. bias
lambdas
## [1] 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
## [18] 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33
## [35] 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
## [52] 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67
## [69] 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84
## [86] 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
ridgeBias <- sapply(lambdas, function(lambda) {
solve(t(learnCovars) %*% learnCovars + lambda * diag(x = 1, J, J)) %*% (t(learnCovars) %*%
learnCovars %*% betas) - betas
})
ridgeBias[1:5, 1:5]
## [,1] [,2] [,3] [,4] [,5]
## [1,] 2.220e-16 0.0043868 0.0086359 0.012756 0.0167533
## [2,] 2.776e-16 0.0000688 0.0002449 0.000516 0.0008713
## [3,] -3.886e-16 0.0018894 0.0037006 0.005440 0.0071122
## [4,] 3.331e-16 0.0005066 0.0009525 0.001346 0.0016949
## [5,] -2.776e-16 -0.0022431 -0.0043303 -0.006278 -0.0080987
# variance
ridgeVar <- sapply(lambdas, function(lambda) {
diag(sigma^2 * solve(t(learnCovars) %*% learnCovars + lambda * diag(x = 1,
J, J)) %*% (t(learnCovars) %*% learnCovars) %*% solve(t(learnCovars) %*%
learnCovars + lambda * diag(x = 1, J, J)))
})
ridgeVar[1:5, 1:5]
## [,1] [,2] [,3] [,4] [,5]
## X 1 0.01616 0.01549 0.01487 0.01429 0.01375
## X 2 0.01781 0.01692 0.01611 0.01536 0.01467
## X 3 0.01459 0.01402 0.01349 0.01299 0.01253
## X 4 0.01728 0.01643 0.01566 0.01496 0.01430
## X 5 0.01961 0.01836 0.01724 0.01624 0.01534
# COMPUTE MSE#
ridgeMSE <- (ridgeBias)^2 + ridgeVar
plot(ridgeMSE)
# best lambda for each Beta
minLambda <- apply(ridgeMSE, 1, which.min)
minLambda
## X 1 X 2 X 3 X 4 X 5 X 6 X 7 X 8 X 9 X 10
## 14 49 55 101 101 101 101 69 31 18
matplot(lambdas, t(ridgeBias), type = "l", lty = 1, col = 1:12, xlab = "Ridge lambda",
ylab = "Bias", main = "", cex.lab = 0.6, cex.axis = 0.6)
legend("top", c(paste0("X", c(1:10))), lty = 3, col = 1:12, cex = 0.5, ncol = 5,
box.col = NA)
matplot(lambdas, t(ridgeMSE), type = "l", lty = 1, col = 1:12, xlab = "Ridge lambda",
ylab = "MSE", main = "", cex.lab = 0.6, cex.axis = 0.6)
legend("top", c(paste0("X", c(1:10))), lty = 3, col = 1:12, cex = 0.5, ncol = 5,
box.col = NA)
# e Lasso Bias, variance MSE
# To simulate these parameters for the Lasso Betas, a MonteCarlo model can
# be made to approximate the true distribution.
newX <- learnCovars %*% betas
MonteCarlo <- function(simNum) {
Y_simNum <- rnorm(n = nLearn, mean = newX, sd = 2)
fit.simNum <- glmnet(x = learnCovars, y = Y_simNum, standardize = F, family = "gaussian",
intercept = F, lambda = lambdasLasso, alpha = alphaLasso)
return(fit.simNum$beta[, c(length(lambdas):1)])
}
resultsBetas <- function(results) {
meanBetas <- NULL
biasBetas <- NULL
varBetas <- NULL
for (i in 1:10) {
betaLassoMC <- sapply(results, function(res) res[i, ])
meanBetas_i <- t(apply(betaLassoMC, 1, mean))
biasBetas_i <- meanBetas_i - betas[i]
varBetas_i <- t(apply(betaLassoMC, 1, var))
meanBetas <- rbind(meanBetas, meanBetas_i)
biasBetas <- rbind(biasBetas, biasBetas_i)
varBetas <- rbind(varBetas, varBetas_i)
}
return(list(meanBetas = meanBetas, biasBetas = biasBetas, varBetas = varBetas))
}
MCResults <- sapply(c(1:1000), MonteCarlo) #create monte carlo model
MCLassoModelBetas <- resultsBetas(MCResults)
str(MCLassoModelBetas)
## List of 3
## $ meanBetas: num [1:10, 1:101] -0.408346 -0.304515 -0.208001 -0.088569 -0.000828 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:101] "s100" "s99" "s98" "s97" ...
## $ biasBetas: num [1:10, 1:101] -0.008346 -0.004515 -0.008001 0.011431 -0.000828 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:101] "s100" "s99" "s98" "s97" ...
## $ varBetas : num [1:10, 1:101] 0.0642 0.0695 0.0559 0.0705 0.0697 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:101] "s100" "s99" "s98" "s97" ...
meanBetas <- MCLassoModelBetas$meanBetas
biasBetas <- MCLassoModelBetas$biasBetas
varBetas <- MCLassoModelBetas$varBetas
mseBetas <- (biasBetas)^2 + varBetas
mseBetas
## s100 s99 s98 s97 s96 s95 s94 s93
## [1,] 0.06431 0.06289 0.06165 0.06058 0.05967 0.05888 0.05819 0.05761
## [2,] 0.06951 0.06680 0.06441 0.06232 0.06041 0.05869 0.05717 0.05585
## [3,] 0.05595 0.05342 0.05108 0.04893 0.04693 0.04506 0.04334 0.04177
## [4,] 0.07058 0.06594 0.06164 0.05767 0.05402 0.05066 0.04758 0.04476
## [5,] 0.06969 0.06333 0.05762 0.05247 0.04783 0.04366 0.03993 0.03656
## [6,] 0.08283 0.07571 0.06922 0.06329 0.05788 0.05293 0.04846 0.04441
## [7,] 0.06644 0.06241 0.05867 0.05519 0.05198 0.04900 0.04626 0.04375
## [8,] 0.07062 0.06665 0.06304 0.05977 0.05682 0.05415 0.05174 0.04958
## [9,] 0.10415 0.09885 0.09409 0.08979 0.08598 0.08256 0.07955 0.07694
## [10,] 0.08837 0.08478 0.08161 0.07879 0.07635 0.07418 0.07234 0.07079
## s92 s91 s90 s89 s88 s87 s86 s85
## [1,] 0.05712 0.05673 0.05640 0.05617 0.05603 0.05595 0.05592 0.05591
## [2,] 0.05467 0.05360 0.05263 0.05179 0.05103 0.05033 0.04969 0.04912
## [3,] 0.04034 0.03903 0.03782 0.03672 0.03571 0.03478 0.03394 0.03319
## [4,] 0.04217 0.03979 0.03759 0.03555 0.03369 0.03197 0.03042 0.02898
## [5,] 0.03355 0.03086 0.02846 0.02632 0.02440 0.02266 0.02111 0.01971
## [6,] 0.04070 0.03732 0.03425 0.03144 0.02888 0.02654 0.02440 0.02247
## [7,] 0.04143 0.03931 0.03735 0.03555 0.03387 0.03232 0.03089 0.02959
## [8,] 0.04763 0.04587 0.04425 0.04276 0.04139 0.04012 0.03896 0.03789
## [9,] 0.07460 0.07253 0.07064 0.06896 0.06743 0.06603 0.06478 0.06369
## [10,] 0.06949 0.06840 0.06744 0.06664 0.06596 0.06534 0.06485 0.06447
## s84 s83 s82 s81 s80 s79 s78 s77
## [1,] 0.05595 0.05602 0.05613 0.05628 0.05650 0.05676 0.05706 0.05739
## [2,] 0.04860 0.04811 0.04766 0.04725 0.04687 0.04653 0.04620 0.04589
## [3,] 0.03250 0.03188 0.03129 0.03076 0.03027 0.02984 0.02946 0.02912
## [4,] 0.02765 0.02640 0.02525 0.02420 0.02320 0.02229 0.02146 0.02068
## [5,] 0.01847 0.01733 0.01631 0.01539 0.01457 0.01381 0.01312 0.01248
## [6,] 0.02072 0.01912 0.01767 0.01635 0.01513 0.01402 0.01300 0.01207
## [7,] 0.02840 0.02729 0.02627 0.02534 0.02449 0.02372 0.02301 0.02235
## [8,] 0.03694 0.03607 0.03529 0.03458 0.03394 0.03337 0.03286 0.03240
## [9,] 0.06272 0.06184 0.06108 0.06041 0.05984 0.05930 0.05883 0.05841
## [10,] 0.06416 0.06393 0.06375 0.06366 0.06366 0.06372 0.06382 0.06395
## s76 s75 s74 s73 s72 s71 s70 s69
## [1,] 0.05776 0.05818 0.058637 0.05913 0.059646 0.060198 0.060790 0.061395
## [2,] 0.04559 0.04531 0.045029 0.04477 0.044524 0.044287 0.044059 0.043834
## [3,] 0.02883 0.02856 0.028320 0.02811 0.027932 0.027771 0.027628 0.027507
## [4,] 0.01997 0.01932 0.018721 0.01817 0.017667 0.017194 0.016752 0.016342
## [5,] 0.01188 0.01132 0.010804 0.01032 0.009873 0.009445 0.009043 0.008668
## [6,] 0.01121 0.01042 0.009689 0.00900 0.008356 0.007754 0.007192 0.006673
## [7,] 0.02174 0.02118 0.020656 0.02016 0.019701 0.019271 0.018866 0.018487
## [8,] 0.03199 0.03163 0.031331 0.03106 0.030823 0.030609 0.030413 0.030234
## [9,] 0.05805 0.05773 0.057453 0.05723 0.057057 0.056908 0.056777 0.056664
## [10,] 0.06410 0.06428 0.064496 0.06475 0.065026 0.065336 0.065679 0.066053
## s68 s67 s66 s65 s64 s63 s62
## [1,] 0.062035 0.062700 0.063394 0.064115 0.064835 0.065565 0.066311
## [2,] 0.043620 0.043417 0.043222 0.043031 0.042833 0.042639 0.042450
## [3,] 0.027406 0.027321 0.027256 0.027206 0.027166 0.027136 0.027119
## [4,] 0.015977 0.015635 0.015313 0.015005 0.014718 0.014447 0.014191
## [5,] 0.008312 0.007982 0.007678 0.007393 0.007124 0.006865 0.006617
## [6,] 0.006195 0.005756 0.005350 0.004970 0.004620 0.004298 0.004007
## [7,] 0.018128 0.017787 0.017468 0.017168 0.016876 0.016594 0.016327
## [8,] 0.030077 0.029933 0.029800 0.029678 0.029565 0.029458 0.029358
## [9,] 0.056581 0.056521 0.056480 0.056458 0.056447 0.056456 0.056479
## [10,] 0.066448 0.066871 0.067311 0.067765 0.068230 0.068731 0.069255
## s61 s60 s59 s58 s57 s56 s55
## [1,] 0.067067 0.067838 0.068633 0.069458 0.070295 0.071143 0.072006
## [2,] 0.042266 0.042091 0.041921 0.041755 0.041591 0.041434 0.041285
## [3,] 0.027122 0.027135 0.027166 0.027198 0.027235 0.027284 0.027341
## [4,] 0.013945 0.013712 0.013491 0.013283 0.013085 0.012893 0.012710
## [5,] 0.006380 0.006153 0.005937 0.005730 0.005532 0.005341 0.005156
## [6,] 0.003734 0.003479 0.003240 0.003017 0.002809 0.002616 0.002435
## [7,] 0.016069 0.015826 0.015596 0.015378 0.015168 0.014965 0.014770
## [8,] 0.029268 0.029189 0.029127 0.029074 0.029033 0.029005 0.028989
## [9,] 0.056508 0.056546 0.056604 0.056685 0.056778 0.056888 0.057024
## [10,] 0.069804 0.070378 0.070976 0.071597 0.072226 0.072877 0.073563
## s54 s53 s52 s51 s50 s49 s48
## [1,] 0.072877 0.073764 0.074662 0.075554 0.076457 0.077380 0.078314
## [2,] 0.041147 0.041015 0.040891 0.040769 0.040654 0.040548 0.040447
## [3,] 0.027408 0.027481 0.027565 0.027656 0.027756 0.027859 0.027972
## [4,] 0.012534 0.012367 0.012209 0.012056 0.011910 0.011775 0.011647
## [5,] 0.004978 0.004806 0.004639 0.004479 0.004322 0.004168 0.004020
## [6,] 0.002266 0.002109 0.001961 0.001824 0.001697 0.001579 0.001469
## [7,] 0.014580 0.014397 0.014218 0.014048 0.013889 0.013736 0.013586
## [8,] 0.028978 0.028968 0.028965 0.028967 0.028972 0.028988 0.029014
## [9,] 0.057179 0.057349 0.057543 0.057755 0.057974 0.058208 0.058455
## [10,] 0.074275 0.075001 0.075739 0.076494 0.077266 0.078065 0.078886
## s47 s46 s45 s44 s43 s42 s41
## [1,] 0.079257 0.080195 0.081138 0.082097 0.083072 0.0840561 0.0850527
## [2,] 0.040349 0.040259 0.040177 0.040109 0.040049 0.0399992 0.0399625
## [3,] 0.028099 0.028236 0.028379 0.028525 0.028675 0.0288249 0.0289747
## [4,] 0.011526 0.011412 0.011306 0.011208 0.011113 0.0110236 0.0109397
## [5,] 0.003879 0.003742 0.003609 0.003479 0.003353 0.0032327 0.0031154
## [6,] 0.001368 0.001273 0.001186 0.001105 0.001029 0.0009594 0.0008957
## [7,] 0.013440 0.013300 0.013164 0.013032 0.012904 0.0127775 0.0126553
## [8,] 0.029052 0.029095 0.029146 0.029202 0.029257 0.0293164 0.0293777
## [9,] 0.058717 0.058988 0.059277 0.059577 0.059902 0.0602518 0.0606160
## [10,] 0.079736 0.080608 0.081498 0.082399 0.083317 0.0842329 0.0851670
## s40 s39 s38 s37 s36 s35 s34
## [1,] 0.0860573 0.0870654 0.0880731 0.089091 0.0901111 0.0911393 0.0921601
## [2,] 0.0399356 0.0399242 0.0399225 0.039932 0.0399567 0.0399937 0.0400370
## [3,] 0.0291295 0.0292896 0.0294513 0.029613 0.0297738 0.0299350 0.0300974
## [4,] 0.0108590 0.0107812 0.0107070 0.010634 0.0105634 0.0104938 0.0104267
## [5,] 0.0030016 0.0028912 0.0027841 0.002679 0.0025772 0.0024783 0.0023827
## [6,] 0.0008371 0.0007827 0.0007321 0.000684 0.0006389 0.0005967 0.0005574
## [7,] 0.0125368 0.0124238 0.0123121 0.012203 0.0120954 0.0119917 0.0118929
## [8,] 0.0294449 0.0295171 0.0295949 0.029675 0.0297514 0.0298289 0.0299129
## [9,] 0.0609881 0.0613838 0.0617811 0.062189 0.0625961 0.0630155 0.0634602
## [10,] 0.0861077 0.0870564 0.0880259 0.089013 0.0900010 0.0909954 0.0919864
## s33 s32 s31 s30 s29 s28
## [1,] 0.0931766 0.0941840 0.0951766 0.0961553 0.0971245 0.0980869
## [2,] 0.0400926 0.0401473 0.0402165 0.0402955 0.0403882 0.0404931
## [3,] 0.0302590 0.0304178 0.0305800 0.0307484 0.0309170 0.0310766
## [4,] 0.0103623 0.0103014 0.0102444 0.0101911 0.0101417 0.0100974
## [5,] 0.0022903 0.0022007 0.0021149 0.0020319 0.0019515 0.0018740
## [6,] 0.0005207 0.0004865 0.0004544 0.0004243 0.0003961 0.0003696
## [7,] 0.0117993 0.0117109 0.0116266 0.0115453 0.0114670 0.0113895
## [8,] 0.0300003 0.0300960 0.0301981 0.0303021 0.0304102 0.0305258
## [9,] 0.0639217 0.0643927 0.0648775 0.0653800 0.0658885 0.0663940
## [10,] 0.0929885 0.0940121 0.0950540 0.0960978 0.0971527 0.0982161
## s27 s26 s25 s24 s23 s22
## [1,] 0.0990512 0.1000130 0.1009719 0.1019312 0.1028864 0.1038519
## [2,] 0.0406137 0.0407467 0.0408867 0.0410345 0.0411932 0.0413631
## [3,] 0.0312338 0.0313908 0.0315506 0.0317109 0.0318752 0.0320414
## [4,] 0.0100557 0.0100176 0.0099800 0.0099444 0.0099084 0.0098747
## [5,] 0.0017991 0.0017272 0.0016582 0.0015914 0.0015272 0.0014652
## [6,] 0.0003447 0.0003214 0.0002997 0.0002794 0.0002603 0.0002425
## [7,] 0.0113154 0.0112444 0.0111763 0.0111105 0.0110458 0.0109838
## [8,] 0.0306452 0.0307719 0.0309042 0.0310363 0.0311707 0.0313087
## [9,] 0.0668945 0.0674014 0.0679074 0.0684161 0.0689188 0.0694171
## [10,] 0.0992904 0.1003572 0.1014383 0.1025208 0.1036114 0.1047093
## s21 s20 s19 s18 s17 s16 s15
## [1,] 0.1048147 0.1057846 0.1067562 0.1077419 0.108730 0.1097109 0.1106971
## [2,] 0.0415480 0.0417413 0.0419503 0.0421687 0.042388 0.0426094 0.0428413
## [3,] 0.0322107 0.0323856 0.0325614 0.0327375 0.032916 0.0330942 0.0332670
## [4,] 0.0098421 0.0098130 0.0097850 0.0097588 0.009735 0.0097153 0.0096969
## [5,] 0.0014055 0.0013480 0.0012930 0.0012403 0.001190 0.0011412 0.0010943
## [6,] 0.0002258 0.0002098 0.0001946 0.0001802 0.000166 0.0001527 0.0001402
## [7,] 0.0109261 0.0108717 0.0108206 0.0107710 0.010721 0.0106710 0.0106237
## [8,] 0.0314467 0.0315878 0.0317275 0.0318645 0.032005 0.0321517 0.0323004
## [9,] 0.0699050 0.0703976 0.0708943 0.0713953 0.071893 0.0723718 0.0728452
## [10,] 0.1058220 0.1069353 0.1080605 0.1091695 0.110253 0.1113498 0.1124381
## s14 s13 s12 s11 s10 s9
## [1,] 0.1116835 0.1126770 0.1136850 1.147e-01 1.157e-01 1.167e-01
## [2,] 0.0430812 0.0433359 0.0435989 4.387e-02 4.415e-02 4.445e-02
## [3,] 0.0334421 0.0336112 0.0337830 3.396e-02 3.413e-02 3.430e-02
## [4,] 0.0096796 0.0096660 0.0096523 9.639e-03 9.626e-03 9.615e-03
## [5,] 0.0010490 0.0010055 0.0009633 9.223e-04 8.818e-04 8.430e-04
## [6,] 0.0001286 0.0001179 0.0001076 9.803e-05 8.912e-05 8.082e-05
## [7,] 0.0105776 0.0105313 0.0104853 1.044e-02 1.040e-02 1.036e-02
## [8,] 0.0324562 0.0326104 0.0327631 3.292e-02 3.307e-02 3.322e-02
## [9,] 0.0733239 0.0737958 0.0742610 7.471e-02 7.517e-02 7.562e-02
## [10,] 0.1135105 0.1145886 0.1156747 1.168e-01 1.179e-01 1.189e-01
## s8 s7 s6 s5 s4 s3
## [1,] 1.176e-01 1.186e-01 1.195e-01 1.204e-01 1.214e-01 0.1222780
## [2,] 4.476e-02 4.509e-02 4.543e-02 4.578e-02 4.614e-02 0.0465072
## [3,] 3.446e-02 3.462e-02 3.478e-02 3.494e-02 3.509e-02 0.0352299
## [4,] 9.606e-03 9.599e-03 9.590e-03 9.582e-03 9.575e-03 0.0095676
## [5,] 8.058e-04 7.696e-04 7.349e-04 7.014e-04 6.683e-04 0.0006363
## [6,] 7.307e-05 6.579e-05 5.896e-05 5.258e-05 4.656e-05 0.0000410
## [7,] 1.032e-02 1.029e-02 1.025e-02 1.022e-02 1.018e-02 0.0101491
## [8,] 3.338e-02 3.354e-02 3.371e-02 3.387e-02 3.403e-02 0.0341827
## [9,] 7.606e-02 7.648e-02 7.690e-02 7.732e-02 7.775e-02 0.0781810
## [10,] 1.200e-01 1.210e-01 1.221e-01 1.231e-01 1.241e-01 0.1251157
## s2 s1 s0
## [1,] 1.232e-01 0.1240407 1.249e-01
## [2,] 4.689e-02 0.0472877 4.769e-02
## [3,] 3.538e-02 0.0355234 3.567e-02
## [4,] 9.561e-03 0.0095572 9.555e-03
## [5,] 6.053e-04 0.0005752 5.463e-04
## [6,] 3.591e-05 0.0000313 2.713e-05
## [7,] 1.011e-02 0.0100806 1.005e-02
## [8,] 3.433e-02 0.0344730 3.461e-02
## [9,] 7.861e-02 0.0790311 7.944e-02
## [10,] 1.261e-01 0.1271145 1.281e-01
lambdaLassoMC <- apply(mseBetas, 1, which.min)
lambdaLassoMC
## [1] 16 63 39 101 101 101 101 49 37 21
for (i in 1:length(lambdaLassoMC)) {
x <- paste("Beta_i=", mseBetas[i, lambdaLassoMC[i]], "maxed by lambda",
lambdaLassoMC[i])
print(x)
}
## [1] "Beta_i= 0.0559114123433861 maxed by lambda 16"
## [1] "Beta_i= 0.0399224617085651 maxed by lambda 63"
## [1] "Beta_i= 0.027118889969228 maxed by lambda 39"
## [1] "Beta_i= 0.00955513332045952 maxed by lambda 101"
## [1] "Beta_i= 0.000546275897290599 maxed by lambda 101"
## [1] "Beta_i= 2.71292588924431e-05 maxed by lambda 101"
## [1] "Beta_i= 0.0100479533432973 maxed by lambda 101"
## [1] "Beta_i= 0.0289647464793083 maxed by lambda 49"
## [1] "Beta_i= 0.0564467569740846 maxed by lambda 37"
## [1] "Beta_i= 0.0636592670456925 maxed by lambda 21"
matplot(lambdas, t(biasBetas), type = "l", lty = 1, col = 1:12, xlab = "Lasso lambda",
ylab = "Bias", main = "", cex.lab = 0.6, cex.axis = 0.6)
legend("top", c(paste("X", c(1:10))), lty = 3, col = 1:12, cex = 0.5, ncol = 5,
box.col = NA)
matplot(lambdas, t(varBetas), type = "l", lty = 1, col = 1:12, xlab = "Lasso lambda",
ylab = "Variance", main = "", cex.lab = 0.6, cex.axis = 0.6)
legend("top", c(paste("X", c(1:10))), lty = 3, col = 1:12, cex = 0.5, ncol = 5,
box.col = NA)
matplot(lambdas, t(mseBetas), type = "l", lty = 1, col = 1:12, xlab = "Lasso lambdas",
ylab = "MSE", main = "", cex.lab = 0.6, cex.axis = 0.6)
legend("top", c(paste0("X", c(1:10))), lty = 3, col = 1:12, cex = 0.5, ncol = 5,
box.col = NA)
