Dades gèniques (403 gens per a 30 mostres) de l’estudi de microarrays de Lu et al. (2004).
lu2004\$`x`: $30\times 403$ matriu que contè el logaritme de les dades. Les files corresponen a les mostres per pacient i les columnes als gens.
lu2004\$y: edat de cada pacient.
The original data are available from the NIH Gene Expression Omnibus and are described in Lu et al. (2004). The selected 403 genes result from prescreening and preprocessing as described in Zuber and Strimmer (2011).
Predir l'edat d'un pacient a través de les variables, que en aquest cas són els gens de les mostres de cervell humà.
Comencem descarregant el paquet care que inclou el conjunt de dades.
#install.packages("care",dependencies=TRUE,repos="https://cloud.r-project.org")
require(care)
data(lu2004)
str(lu2004)
List of 2 $ x: num [1:30, 1:403] 6.88 7.14 6.3 6.61 6.64 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : chr [1:30] "GSM27015" "GSM27016" "GSM27017" "GSM27018" ... .. ..$ : chr [1:403] "1007_s_at" "1009_at" "1040_s_at" "1074_at" ... $ y: num [1:30] 26 26 27 29 30 36 37 38 40 42 ...
head(lu2004)
| 1007_s_at | 1009_at | 1040_s_at | 1074_at | 1090_f_at | 1130_at | 113_i_at | 1158_s_at | 1161_at | 1163_at | ... | 942_at | 953_g_at | 955_at | 959_at | 976_s_at | 996_at | AFFX-HSAC07/X00351_5_at | AFFX-HSAC07/X00351_M_at | AFFX-HUMGAPDH/M33197_5_at | AFFX-HUMISGF3A/M97935_5_at | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| GSM27015 | 6.876809 | 7.982007 | 5.477952 | 6.313889 | 6.326945 | 6.998656 | 4.820709 | 8.187118 | 7.994106 | 4.830272 | ... | 5.642811 | 6.212414 | 7.980304 | 5.597002 | 4.936364 | 6.247186 | 8.533037 | 8.674293 | 8.753856 | 5.531732 |
| GSM27016 | 7.142606 | 8.011060 | 5.413377 | 6.051524 | 6.100904 | 7.230498 | 5.485154 | 8.279875 | 8.296300 | 4.012053 | ... | 5.284350 | 6.548317 | 8.048235 | 5.551299 | 4.977664 | 6.452093 | 8.631309 | 8.698900 | 8.631343 | 3.522273 |
| GSM27017 | 6.299614 | 7.840191 | 5.418564 | 6.228191 | 6.203052 | 6.588474 | 5.471808 | 7.370445 | 8.396026 | 4.954227 | ... | 4.429384 | 6.102659 | 7.704393 | 5.638212 | 4.840630 | 5.075018 | 8.399171 | 8.497056 | 8.816199 | 4.127104 |
| GSM27018 | 6.611963 | 7.965969 | 5.776317 | 6.150937 | 6.106631 | 6.979006 | 4.870768 | 8.257948 | 8.146726 | 5.040298 | ... | 5.129786 | 6.353547 | 7.972349 | 5.155267 | 5.150988 | 5.174951 | 8.524826 | 8.678832 | 8.633710 | 4.042114 |
| GSM27019 | 6.636269 | 7.914380 | 5.659092 | 5.892539 | 6.071977 | 6.929752 | 5.518220 | 8.181841 | 8.345850 | 5.088466 | ... | 4.561765 | 6.532718 | 7.928276 | 5.399781 | 5.183854 | 5.195171 | 8.686615 | 8.727050 | 8.761970 | 4.140913 |
| GSM27020 | 6.869010 | 7.813583 | 5.258916 | 5.670536 | 6.286869 | 6.969386 | 4.841119 | 8.126205 | 8.163049 | 4.662099 | ... | 4.822980 | 6.205765 | 7.917759 | 4.460732 | 4.423890 | 6.374539 | 8.448844 | 8.581228 | 8.542085 | 3.695177 |
| GSM27021 | 6.616485 | 7.859228 | 5.357845 | 6.123738 | 6.372847 | 6.542662 | 4.567167 | 7.956589 | 7.950185 | 4.991520 | ... | 5.339267 | 5.725417 | 7.877348 | 4.773621 | 4.564296 | 6.315879 | 8.207263 | 8.459058 | 8.597812 | 4.095521 |
| GSM27022 | 6.639573 | 7.967668 | 5.352408 | 5.893609 | 6.384028 | 6.879522 | 4.518702 | 8.188664 | 8.165858 | 4.727874 | ... | 5.424086 | 6.093437 | 7.931795 | 5.320597 | 4.407340 | 5.598226 | 8.299395 | 8.496172 | 8.525705 | 3.755308 |
| GSM27023 | 6.724338 | 7.805308 | 4.970986 | 5.327978 | 6.228870 | 6.447911 | 4.291531 | 7.967044 | 7.909941 | 4.345730 | ... | 5.248055 | 5.765868 | 7.570943 | 4.927725 | 4.035053 | 5.947191 | 8.241832 | 8.423445 | 8.350021 | 3.154495 |
| GSM27024 | 6.478752 | 7.921695 | 5.191901 | 5.553281 | 6.154794 | 6.748245 | 5.634286 | 8.104395 | 8.036441 | 4.207197 | ... | 4.773300 | 5.377466 | 7.752627 | 5.044876 | 4.081153 | 5.805933 | 7.934456 | 8.239741 | 8.496290 | 3.567737 |
| GSM27025 | 7.177798 | 7.842067 | 5.138495 | 5.568356 | 6.236507 | 6.696230 | 4.070192 | 7.684126 | 7.918643 | 4.365450 | ... | 5.570388 | 5.145090 | 7.842165 | 4.973971 | 4.069734 | 6.625033 | 7.659520 | 8.187730 | 8.292954 | 3.539079 |
| GSM27026 | 6.652482 | 7.900714 | 5.833892 | 5.905425 | 6.406875 | 7.276466 | 4.693749 | 7.756238 | 8.077081 | 4.662609 | ... | 4.747459 | 5.663217 | 7.697530 | 5.195637 | 4.904282 | 6.149801 | 7.938146 | 8.211010 | 8.264641 | 3.955194 |
| GSM27027 | 7.032139 | 7.975042 | 4.791650 | 5.808067 | 6.375345 | 7.048213 | 4.636737 | 7.981354 | 8.255540 | 4.330995 | ... | 5.407481 | 5.638917 | 7.976633 | 5.170995 | 3.567466 | 6.320989 | 8.013558 | 8.381532 | 8.421757 | 2.895475 |
| GSM27028 | 6.888724 | 7.780023 | 5.499244 | 5.987534 | 6.418884 | 6.910771 | 5.215691 | 7.981426 | 7.970699 | 4.739019 | ... | 5.565711 | 6.001447 | 7.905685 | 5.257204 | 4.843045 | 6.547021 | 8.222511 | 8.341937 | 8.505847 | 3.637460 |
| GSM27029 | 6.973356 | 7.829380 | 5.293229 | 6.164226 | 6.335795 | 6.688523 | 4.405527 | 7.869279 | 8.122718 | 4.889236 | ... | 5.716021 | 5.781586 | 7.716768 | 5.085489 | 4.867473 | 6.735172 | 8.311347 | 8.441758 | 8.536576 | 3.683249 |
| GSM27030 | 6.638881 | 7.926700 | 5.089009 | 5.768790 | 6.330284 | 6.749905 | 4.935883 | 8.188081 | 8.199843 | 4.354997 | ... | 4.826897 | 5.906138 | 7.604730 | 5.321814 | 3.894922 | 5.786499 | 8.233503 | 8.483028 | 8.400500 | 3.619666 |
| GSM27031 | 7.327426 | 7.972190 | 4.983524 | 5.966772 | 6.332496 | 7.128103 | 4.593984 | 7.928896 | 8.149229 | 4.496569 | ... | 5.361499 | 5.805129 | 7.983911 | 5.456457 | 4.287539 | 7.303150 | 8.401699 | 8.580829 | 8.603841 | 3.622026 |
| GSM27032 | 6.984957 | 7.520430 | 4.812079 | 5.320680 | 6.362221 | 6.260485 | 4.065827 | 7.831284 | 7.716666 | 4.107909 | ... | 5.379856 | 4.693044 | 7.537776 | 4.235656 | 3.858510 | 6.538902 | 7.849035 | 8.271418 | 8.045431 | 3.699119 |
| GSM27033 | 7.053586 | 7.765726 | 5.345335 | 5.869893 | 6.376756 | 6.868551 | 4.488961 | 8.007567 | 7.987840 | 4.518110 | ... | 5.702058 | 5.787958 | 7.851097 | 5.099049 | 4.592603 | 6.685052 | 8.184048 | 8.393252 | 8.522165 | 3.822155 |
| GSM27034 | 7.110712 | 7.785072 | 4.754977 | 5.040621 | 6.384096 | 6.412153 | 4.039004 | 7.414958 | 7.711003 | 3.877846 | ... | 5.700654 | 4.317909 | 7.865898 | 4.260962 | 3.731292 | 6.539303 | 7.528364 | 8.179242 | 7.951785 | 3.265778 |
| GSM27035 | 7.110321 | 7.652693 | 5.119884 | 5.703536 | 6.424190 | 6.549343 | 4.452423 | 7.856653 | 8.062578 | 4.386540 | ... | 5.710735 | 5.030196 | 7.843170 | 4.528162 | 4.149044 | 6.644074 | 7.772285 | 8.085157 | 8.213785 | 3.463674 |
| GSM27036 | 7.205346 | 7.572724 | 5.080889 | 5.066076 | 6.452090 | 5.871946 | 4.501778 | 7.540791 | 7.774284 | 4.030226 | ... | 5.429850 | 5.292390 | 7.538234 | 3.132900 | 4.277017 | 6.510115 | 8.124284 | 8.296147 | 8.465761 | 3.613668 |
| GSM27037 | 7.022476 | 7.906768 | 5.048117 | 6.209937 | 6.293157 | 7.220900 | 4.201443 | 8.042050 | 8.221882 | 4.412128 | ... | 5.737936 | 5.273271 | 7.846081 | 5.466160 | 4.113235 | 6.630718 | 7.991822 | 8.305546 | 8.370860 | 3.391433 |
| GSM27038 | 7.196087 | 7.812872 | 4.823896 | 5.707877 | 6.207690 | 6.877099 | 3.785035 | 7.307712 | 7.886739 | 3.989948 | ... | 6.059315 | 4.799733 | 7.740225 | 5.451858 | 3.806045 | 7.271447 | 7.641334 | 8.082282 | 8.163235 | 2.525953 |
| GSM27039 | 6.941770 | 7.913082 | 5.021766 | 5.932627 | 6.215783 | 6.783972 | 3.771464 | 8.053725 | 7.996741 | 4.635864 | ... | 5.663922 | 5.545275 | 7.853294 | 5.401911 | 4.135370 | 6.426501 | 8.203726 | 8.522348 | 8.523911 | 2.982596 |
| GSM27040 | 6.987601 | 7.899803 | 5.082757 | 5.474726 | 6.342678 | 6.781468 | 4.129400 | 8.098622 | 8.140785 | 4.178244 | ... | 5.509967 | 5.298102 | 7.790679 | 4.817843 | 3.946279 | 6.779529 | 8.052322 | 8.424486 | 8.236995 | 3.294684 |
| GSM27041 | 7.022182 | 7.383001 | 4.730251 | 5.016332 | 6.611315 | 6.180188 | 4.716381 | 7.477485 | 7.563154 | 4.150983 | ... | 5.691197 | 5.532666 | 7.684844 | 4.099669 | 4.388507 | 6.815761 | 8.081481 | 8.177564 | 8.340007 | 2.855027 |
| GSM27042 | 7.085332 | 7.652114 | 5.136040 | 5.449926 | 6.368199 | 6.168445 | 4.458586 | 7.676163 | 7.904667 | 4.359543 | ... | 5.373541 | 5.509761 | 7.552573 | 4.646149 | 4.410313 | 5.939537 | 8.426382 | 8.572083 | 8.518986 | 3.801822 |
| GSM27043 | 7.082960 | 7.257334 | 4.730092 | 4.567300 | 6.888915 | 5.586878 | 4.513064 | 6.352950 | 7.420202 | 3.992913 | ... | 5.212656 | 4.681557 | 6.909074 | 3.936920 | 4.162315 | 7.045872 | 7.803231 | 8.030257 | 8.459903 | 3.325471 |
| GSM27044 | 7.036588 | 7.356427 | 4.756826 | 5.357609 | 6.595670 | 5.574691 | 3.731132 | 6.691723 | 7.317305 | 3.937137 | ... | 5.906274 | 4.322992 | 7.143570 | 2.372354 | 3.962309 | 6.773539 | 7.245070 | 7.642275 | 8.202115 | 3.338566 |
Abans de tot, hem de veure si falten valors d'alguna variable per tal d'ometre aquestes dades buides.
Utilitzarem la funció is.na() que indica les observacions buides i la funció sum() que ens dona el nombre d'elements que falten.
x<-lu2004$x
y<-as.numeric(lu2004$y)
dim(x)
sum(is.na(lu2004$y))
A continuació apliquem el mètode de regressió lineal de Mínims Quadrats Ordinaris (MQO), on y serà la variable aleatòria resposta, i totes les altres formaran la matriu $X$ de predictors.
Q<- t(x) %*% x #t() transposada
det(Q) #Calcuar determinant
Q1<-solve(Q) #solve() inversa
round(Q1,3) #Printar i aproximar matriu Q^-1, 3 decimals
Error in solve.default(Q): sistema es computacionalmente singular: número de condición recíproco = 4.68919e-21 Traceback: 1. solve(Q) 2. solve.default(Q)
No es pot aplicar MQO ja que la matriu $Q=X^{t}X$ és singular ($N<p$) i per tant no es poden calcular els estimadors.
round(kappa(x),2)
Com podem veure, el número de condició és molt elevat.
Com que la restricció del ridge/lasso/Elastic Net tracta tots els coeficients per igual, normalment té sentit que tots els elements de X estiguin en les mateixes unitats. Si no, sense cap pèrdua de generalitat, centralitzarem i estandarditzarem els predictors per tal que cadascuna de les columnes de X tingui mitjana 0 i variància 1.
xm<-apply(x,2,mean)
xc<-sweep(x,2,xm,"-")
xs<-apply(xc,2,sd)
x0<-sweep(xc,2,xs,"/")
str(x0)
round(max(abs(apply(x0,2,mean))),6)
round(max(abs(apply(x0,2,var))),6)
num [1:30, 1:403] -0.152 0.927 -2.493 -1.226 -1.127 ... - attr(*, "dimnames")=List of 2 ..$ : chr [1:30] "GSM27015" "GSM27016" "GSM27017" "GSM27018" ... ..$ : chr [1:403] "1007_s_at" "1009_at" "1040_s_at" "1074_at" ...
Arribem a la conclusió que com hi ha multicol·linealitat entre les variables del model no podem aplicar MQO. Per tant, utilitzarem la regressió Ridge per estimar els coeficients del model.
Preparem les dades en forma de matriu per poder utilitzar la funció glmnet, un paquet que inclou les regressions Ridge, Lasso i Elastic Net que farem servir.
x<-x0
#install.packages("glmnet",dependencies=TRUE,repos="https://cloud.r-project.org")
require(glmnet)
Primer definim un ventall de possible valors de $\lambda$:
grid<-exp(seq(10,-3,length=100))
Per obtenir el camí de solucions dels coeficients en funció de $\lambda$ apliquem la funció glmnet() a les dades utilitzant que:
alpha=0 correspon a la regressió Ridge
alpha=1 correspon a la regressió Lasso
lu2004.ridge.01<-glmnet(x,y,alpha=0,lambda=grid)
Ara utilitzem la funció plotper realitzar el gràfic que ens permet veure l'estimació dels coeficientes del model ajustat respecte a:
options(repr.plot.width=20,repr.plot.height=7)
old.par<-par(mfrow=c(1,3))
plot(lu2004.ridge.01,xvar="norm")
abline(h=0,lty=2,col='black')
title('Ridge Regression\n')
plot(lu2004.ridge.01,xvar="lambda")
abline(h=0,lty=2,col='black')
title('Ridge Regression\n')
plot(lu2004.ridge.01,xvar="dev")
abline(h=0,lty=2,col='black')
title('Ridge Regression\n')
Per a calcular el valor de $\lambda$ òptim, primer hem d'establir una "llavor aleatòria" perquè els resultats obtinguts puguin ser reproduïbles.
Després apliquem la funció cv.glmnet() per obtenir el valor de $\lambda$ que minimitza l'error de predicció, mitjançant la validació encreuada.
set.seed(24025)
lu2004.ridge.01.cv<-cv.glmnet(x,y,alpha=0,lambda=grid)
Tornem a utilitzar la funció plot per veure gràficament els resultats de la validació encreuada. Podem veure dues marques:
bestlam.ridge1<-lu2004.ridge.01.cv$lambda.min
lam1se.ridge1<-lu2004.ridge.01.cv$lambda.1se
options(repr.plot.width=7, repr.plot.height=7)
plot(lu2004.ridge.01.cv)
abline(v=log(bestlam.ridge1),lwd=3,col="cyan")
round(log(bestlam.ridge1),2)
round(log(lam1se.ridge1),2)
Finalment utilitzem la funció predict() per reajustar els coeficients utilitzant el valor de la $\lambda$ resultant de la validació encreuada.
# We extract now the regression coefficients with the 'predict' function
coeff.01.ridge<-as.numeric(predict(lu2004.ridge.01,s=bestlam.ridge1,type="coefficients"))
round(coeff.01.ridge,2)
Com podem veure, els coeficients han quedat restringits cap a zero i a canvi d'afegir biaix hem reduït la variància. Tanmateix cap dels coeficients és igual a zero, ja que la regressió Ridge no fa selecció de variables.
Repetim el mateix procediment, però en aquest cas, estimant el model mitjançant la regressió Lasso.
Primer definim un ventall de possible valors de $\lambda$:
grid<-exp(seq(10,-5,length=100))
#str(grid)
Per obtenir el camí de solucions dels coeficients en funció de $\lambda$ apliquem la funció glmnet() a les dades utilitzant que alpha=1 correspon a la regressió Lasso.
lu2004.lasso.01<-glmnet(x,y,alpha=1,lambda=grid)
Ara utilitzem la funció plotper realitzar el gràfic que ens permet veure l'estimació dels coeficientes del model ajustat respecte a:
options(repr.plot.width=20,repr.plot.height=7)
old.par<-par(mfrow=c(1,3))
plot(lu2004.lasso.01,xvar="norm")
abline(h=0,lty=2,col='black')
title('Lasso Regression\n')
plot(lu2004.lasso.01,xvar="lambda")
abline(h=0,lty=2,col='black')
title('Lasso Regression\n')
plot(lu2004.lasso.01,xvar="dev")
abline(h=0,lty=2,col='black')
title('Lasso Regression\n')
Per a calcular el valor de $\lambda$ òptim, tornem a establir la mateixa "llavor aleatòria" perquè els resultats obtinguts puguin ser reproduïbles.
Després apliquem la funció cv.glmnet() per obtenir el valor de $\lambda$ que minimitza l'error de predicció, mitjançant la validació encreuada.
set.seed(24025)
lu2004.lasso.01.cv<-cv.glmnet(x,y,alpha=1,lambda=grid)
Tornem a utilitzar la funció plot per veure gràficament els resultats de la validació encreuada. Podem veure dues marques:
bestlam.lasso1<-lu2004.lasso.01.cv$lambda.min
lam1se.lasso1<-lu2004.lasso.01.cv$lambda.1se
options(repr.plot.width=7, repr.plot.height=7)
plot(lu2004.lasso.01.cv)
abline(v=log(bestlam.lasso1),lwd=3,col="cyan")
round(log(bestlam.lasso1),2)
round(log(lam1se.lasso1),2)
Finalment utilitzem la funció predict() per reajustar els coeficients utilitzant el valor de la $\lambda$ resultant de la validació encreuada.
# We extract now the regression coefficients with the 'predict' function
coeff.01.lasso<-as.numeric(predict(lu2004.lasso.01,s=bestlam.lasso1,type="coefficients"))
lu2004.Intercept<-coeff.01.lasso[1]
round(coeff.01.lasso,2)
Idx.Non.null.coeff<-which(coeff.01.lasso!=0)
Idx.Non.null.coeff
length(Idx.Non.null.coeff)
lu2004.vars<-dimnames(lu2004$x)[[2]]
str(lu2004.vars)
chr [1:403] "1007_s_at" "1009_at" "1040_s_at" "1074_at" "1090_f_at" ...
lu2004.result<-rbind(c("Int.",lu2004.vars[Idx.Non.null.coeff]),
round(c(lu2004.Intercept,coeff.01.lasso[Idx.Non.null.coeff]),2))
dim(lu2004.result)<-c(8,9)
lu2004.result
| Int. | 296_at | 33299_at | 34887_at | 36273_at | 37812_at | 39395_at | 40608_at | 41494_at |
| 60.27 | -0.92 | -0.05 | 1.76 | -0.12 | -1.49 | 2.55 | 2.42 | 1.32 |
| 1007_s_at | 317_at | 33508_at | 353_at | 36823_at | 37908_at | 39850_at | 41086_at | 871_s_at |
| 60.27 | -5.25 | -4.62 | -0.03 | -0.21 | -3.2 | -0.32 | -0.31 | 3.38 |
| 1593_at | 32215_i_at | 33742_f_at | 35494_at | 37301_at | 38515_at | 39983_at | 41479_s_at | 959_at |
| 1.8 | 0 | -1.76 | 1.41 | -0.27 | 3.6 | -1.44 | -0.03 | 0 |
| 2083_at | 32916_at | 34272_at | 35668_at | 37712_g_at | 38627_at | 40098_at | 41483_s_at | NA |
| 0.01 | 5.35 | -1.47 | -4.19 | -0.08 | -1.52 | 0.08 | -0.49 | -2.32 |
Veiem que només 35 dels 403 coeficients de la regressió són diferents de zero, és a dir, les altres variables associades a coeficients iguals a zero no estaran en el model.
Hem comprobat que el lasso realitza selecció de variables.
El procediment Elastic Net pren valors de alpha entre $0$ i $1$. Anem a veure com canvia el model per a alphes diferents.
Definim un ventall de possible valors de $\lambda$:
grid<-exp(seq(9,-20,length=100))
Per obtenir el camí de solucions dels coeficients en funció de $\lambda$ apliquem la funció glmnet() a les dades utilitzant:
alpha=0.001
alpha=0.3
alpha=0.5
alpha=0.999
lu2004.net.01<-glmnet(x,y,alpha=0.001,lambda=grid)
lu2004.net.02<-glmnet(x,y,alpha=0.3,lambda=grid)
lu2004.net.03<-glmnet(x,y,alpha=0.5,lambda=grid)
lu2004.net.04<-glmnet(x,y,alpha=0.999,lambda=grid)
Ara utilitzem la funció plotper realitzar el gràfic que ens permet veure l'estimació dels coeficientes del model ajustat per a les diferents alpha:
options(repr.plot.width=20,repr.plot.height=7)
old.par<-par(mfrow=c(1,4))
plot(lu2004.net.01,xvar="norm")
abline(h=0,lty=2,col='black')
title('alpha=0.001\n')
plot(lu2004.net.02,xvar="norm")
abline(h=0,lty=2,col='black')
title('alpha=0.3\n')
plot(lu2004.net.03,xvar="norm")
abline(h=0,lty=2,col='black')
title('alpha=0.5\n')
plot(lu2004.net.04,xvar="norm")
abline(h=0,lty=2,col='black')
title('alpha=0.999\n')
Com podem observar, la primera gràfica (alpha=0.001) és molt semblant al camí de les solucions que proporciona la regressió Ridge, ja que el coeficient alpha està molt a prop de 0. Per altra banda, l'última gràfica (alpha=0.999) és quasi idèntica al camí de les solucions que ens dona la regressió Lasso, perquè el coeficient alpha és quasi 1.
Tornem a establir la mateixa "llavor aleatòria" per a que els resultats obtinguts siguin reproduïbles.
Després apliquem la funció cv.glmnet() per obtenir el valor de $\lambda$ que minimitza l'error de predicció, mitjançant la validació encreuada, aplicada a cada model.
set.seed(23025)
lu2004.net.01.cv<-cv.glmnet(x,y,alpha=0.001,lambda=grid)
lu2004.net.02.cv<-cv.glmnet(x,y,alpha=0.3,lambda=grid)
lu2004.net.03.cv<-cv.glmnet(x,y,alpha=0.5,lambda=grid)
lu2004.net.04.cv<-cv.glmnet(x,y,alpha=0.999,lambda=grid)
Tornem a utilitzar la funció plot per veure gràficament els resultats de la validació encreuada. Podem veure dues marques:
options(repr.plot.width=15,repr.plot.height=15)
old.par<-par(mfrow=c(2,2))
bestlam.net1<-lu2004.net.01.cv$lambda.min
lam1se.net1<-lu2004.net.01.cv$lambda.1se
plot(lu2004.net.01.cv)
abline(v=log(bestlam.net1),lwd=3,col="cyan")
round(log(bestlam.net1),2)
round(log(lam1se.net1),2)
bestlam.net2<-lu2004.net.02.cv$lambda.min
lam1se.net2<-lu2004.net.02.cv$lambda.1se
plot(lu2004.net.02.cv)
abline(v=log(bestlam.net2),lwd=3,col="cyan")
round(log(bestlam.net2),2)
round(log(lam1se.net2),2)
bestlam.net3<-lu2004.net.03.cv$lambda.min
lam1se.net3<-lu2004.net.03.cv$lambda.1se
plot(lu2004.net.03.cv)
abline(v=log(bestlam.net3),lwd=3,col="cyan")
round(log(bestlam.net3),2)
round(log(lam1se.net3),2)
bestlam.net4<-lu2004.net.04.cv$lambda.min
lam1se.net4<-lu2004.net.04.cv$lambda.1se
plot(lu2004.net.04.cv)
abline(v=log(bestlam.net4),lwd=3,col="cyan")
round(log(bestlam.net4),2)
round(log(lam1se.net4),2)
Finalment utilitzem la funció predict() per reajustar els coeficients utilitzant el valor de la $\lambda$ resultant de la validació encreuada.
coeff.01.net<-as.numeric(predict(lu2004.net.01,s=bestlam.net1,type="coefficients"))
coeff.02.net<-as.numeric(predict(lu2004.net.02,s=bestlam.net2,type="coefficients"))
coeff.03.net<-as.numeric(predict(lu2004.net.03,s=bestlam.net3,type="coefficients"))
coeff.04.net<-as.numeric(predict(lu2004.net.04,s=bestlam.net4,type="coefficients"))
Idx.Non.null.coeff1<-sum(coeff.01.net!=0)
Idx.Non.null.coeff1
Idx.Non.null.coeff2<-sum(coeff.02.net!=0)
Idx.Non.null.coeff2
Idx.Non.null.coeff3<-sum(coeff.03.net!=0)
Idx.Non.null.coeff3
Idx.Non.null.coeff4<-sum(coeff.04.net!=0)
Idx.Non.null.coeff4
Similars als coeficients de la regressió lasso, amb alguns coeficients nuls, que resolen les limitacions del lasso.
A continuació podem veure els diferents gràfics que ens permeten veure l'estimació dels coeficients del model ajustat pels diferents alpha, incloent-hi les regressions Ridge i Lasso.
options(repr.plot.width=14,repr.plot.height=10)
old.par<-par(mfrow=c(2,3))
plot(lu2004.ridge.01)
abline(h=0,lty=2,col='black')
title('alpha=0 (Ridge)\n')
plot(lu2004.net.01,xvar="norm")
abline(h=0,lty=2,col='black')
title('alpha=0.001\n')
plot(lu2004.net.02,xvar="norm")
abline(h=0,lty=2,col='black')
title('alpha=0.3\n')
plot(lu2004.net.03,xvar="norm")
abline(h=0,lty=2,col='black')
title('alpha=0.5\n')
plot(lu2004.net.04,xvar="norm")
abline(h=0,lty=2,col='black')
title('alpha=0.999\n')
plot(lu2004.lasso.01)
abline(h=0,lty=2,col='black')
title('alpha=1 (Lasso)\n')
Finalment, podem comparar els coeficients obtinguts mitjançant la lambda que minimitza la mitjana de l'error per a la regressió Lasso.
print('Coeficients regressió Ridge')
coeff.01.ridge<-as.numeric(predict(lu2004.ridge.01,s=bestlam.lasso1,type="coefficients"))
Idx.Non.null.coeff0<-sum(coeff.01.ridge!=0)
Idx.Non.null.coeff0
print('Coeficients regressió Elastic Net')
coeff.01.net<-as.numeric(predict(lu2004.net.01,s=bestlam.lasso1,type="coefficients"))
coeff.02.net<-as.numeric(predict(lu2004.net.02,s=bestlam.lasso1,type="coefficients"))
coeff.03.net<-as.numeric(predict(lu2004.net.03,s=bestlam.lasso1,type="coefficients"))
coeff.04.net<-as.numeric(predict(lu2004.net.04,s=bestlam.lasso1,type="coefficients"))
Idx.Non.null.coeff1<-sum(coeff.01.net!=0)
Idx.Non.null.coeff1
Idx.Non.null.coeff2<-sum(coeff.02.net!=0)
Idx.Non.null.coeff2
Idx.Non.null.coeff3<-sum(coeff.03.net!=0)
Idx.Non.null.coeff3
Idx.Non.null.coeff4<-sum(coeff.04.net!=0)
Idx.Non.null.coeff4
print('Coeficients regressió Lasso')
Idx.Non.null.coeff5<-sum(coeff.01.lasso!=0)
Idx.Non.null.coeff5
[1] "Coeficients regressió Ridge"
[1] "Coeficients regressió Elastic Net"
[1] "Coeficients regressió Lasso"
Veiem que els coeficients provenents del ridge no presenten cap estimació nul·la, mentre que els coeficients de l'elastic net a mesura que s'aproximen a la regressió lasso, començen a tenir alguns coeficients iguals a zero.