Aquestes dades provenen d'un estudi que examina la correlació que hi ha entre el nivell d'un antigen específic de la pròstata i un nombre de mesures clíniques en homes que estaven a punt de rebre una prostatectomia radical.
És un dataframe de 97 files i 9 columnes.
Aquestes dades es componen de la següent informació:
lcavol: log(volum de càncer)lweight: log(pes de la pròstata)age: edadlbph: log(quantitat d'hiperplàsia prostàtica benigna)svi: invasió de vesícules seminalslcp: log(penetració capsular)gleason: puntuació de Gleasonpgg45: percentatge de puntuacions de Gleason 4 o 5lpsa: log(antigen específic de la pròstata)Stamey, T.A., Kabalin, J.N., McNeal, J.E., Johnstone, I.M., Freiha, F., Redwine, E.A. and Yang, N. (1989), "Prostate specific antigen in the diagnosis and treatment of adenocarcinoma of the prostate: II. radical prostatectomy treated patients", Journal of Urology 141(5), 1076–1083.
Predir la variable resposta lpsa, mitjançant totes les altres variables de la base de dades que actuaran com a predictors.
Comencem descarregant el paquet lasso2 que inclou el conjunt de dades Prostate.
#install.packages("lasso2",dependencies=TRUE,repos="https://cloud.r-project.org")
require(lasso2)
data(Prostate)
str(Prostate)
'data.frame': 97 obs. of 9 variables: $ lcavol : num -0.58 -0.994 -0.511 -1.204 0.751 ... $ lweight: num 2.77 3.32 2.69 3.28 3.43 ... $ age : num 50 58 74 58 62 50 64 58 47 63 ... $ lbph : num -1.39 -1.39 -1.39 -1.39 -1.39 ... $ svi : num 0 0 0 0 0 0 0 0 0 0 ... $ lcp : num -1.39 -1.39 -1.39 -1.39 -1.39 ... $ gleason: num 6 6 7 6 6 6 6 6 6 6 ... $ pgg45 : num 0 0 20 0 0 0 0 0 0 0 ... $ lpsa : num -0.431 -0.163 -0.163 -0.163 0.372 ...
head(Prostate)
| lcavol | lweight | age | lbph | svi | lcp | gleason | pgg45 | lpsa | |
|---|---|---|---|---|---|---|---|---|---|
| <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | |
| 1 | -0.5798185 | 2.769459 | 50 | -1.386294 | 0 | -1.386294 | 6 | 0 | -0.4307829 |
| 2 | -0.9942523 | 3.319626 | 58 | -1.386294 | 0 | -1.386294 | 6 | 0 | -0.1625189 |
| 3 | -0.5108256 | 2.691243 | 74 | -1.386294 | 0 | -1.386294 | 7 | 20 | -0.1625189 |
| 4 | -1.2039728 | 3.282789 | 58 | -1.386294 | 0 | -1.386294 | 6 | 0 | -0.1625189 |
| 5 | 0.7514161 | 3.432373 | 62 | -1.386294 | 0 | -1.386294 | 6 | 0 | 0.3715636 |
| 6 | -1.0498221 | 3.228826 | 50 | -1.386294 | 0 | -1.386294 | 6 | 0 | 0.7654678 |
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.
dim(Prostate)
sum(is.na(Prostate$lpsa))
A continuació apliquem el mètode de regressió lineal de Mínims Quadrats Ordinaris (MQO), on lpsa serà la variable aleatòria resposta, i totes les altres formaran la matriu X de predictors.
Prostate.lm.1<-lm(lpsa~., data=Prostate)
summary(Prostate.lm.1)
Call:
lm(formula = lpsa ~ ., data = Prostate)
Residuals:
Min 1Q Median 3Q Max
-1.73316 -0.37133 -0.01702 0.41414 1.63811
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.669399 1.296381 0.516 0.60690
lcavol 0.587023 0.087920 6.677 2.11e-09 ***
lweight 0.454461 0.170012 2.673 0.00896 **
age -0.019637 0.011173 -1.758 0.08229 .
lbph 0.107054 0.058449 1.832 0.07040 .
svi 0.766156 0.244309 3.136 0.00233 **
lcp -0.105474 0.091013 -1.159 0.24964
gleason 0.045136 0.157464 0.287 0.77506
pgg45 0.004525 0.004421 1.024 0.30885
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7084 on 88 degrees of freedom
Multiple R-squared: 0.6548, Adjusted R-squared: 0.6234
F-statistic: 20.86 on 8 and 88 DF, p-value: < 2.2e-16
Abans de donar aquest model per bo, ens hem d'assegurar que els predictors no presenten un problema de multicol·linealitat. Primer vegem la matriu de correlacions en forma numèrica i en forma gràfica:
x<-model.matrix(lpsa~.,Prostate)[,-1]
round(cor(x),2)->mat_cor
mat_cor
| lcavol | lweight | age | lbph | svi | lcp | gleason | pgg45 | |
|---|---|---|---|---|---|---|---|---|
| lcavol | 1.00 | 0.19 | 0.22 | 0.03 | 0.54 | 0.68 | 0.43 | 0.43 |
| lweight | 0.19 | 1.00 | 0.31 | 0.43 | 0.11 | 0.10 | 0.00 | 0.05 |
| age | 0.22 | 0.31 | 1.00 | 0.35 | 0.12 | 0.13 | 0.27 | 0.28 |
| lbph | 0.03 | 0.43 | 0.35 | 1.00 | -0.09 | -0.01 | 0.08 | 0.08 |
| svi | 0.54 | 0.11 | 0.12 | -0.09 | 1.00 | 0.67 | 0.32 | 0.46 |
| lcp | 0.68 | 0.10 | 0.13 | -0.01 | 0.67 | 1.00 | 0.51 | 0.63 |
| gleason | 0.43 | 0.00 | 0.27 | 0.08 | 0.32 | 0.51 | 1.00 | 0.75 |
| pgg45 | 0.43 | 0.05 | 0.28 | 0.08 | 0.46 | 0.63 | 0.75 | 1.00 |
#install.packages("corrplot",dependencies=TRUE,repos="https://cloud.r-project.org")
require(corrplot)
options(repr.plot.width=8,repr.plot.height=8)
corrplot(mat_cor, type="upper", order="hclust", tl.col="black", tl.srt=45)
Veiem que alguns coeficients de correlació són bastant propers a 1, el que indica que hi ha relacions entre els predictors.
Calcularem el FIV, funció que es troba al paquet car:
#install.packages("car",dependencies=TRUE,repos="https://cloud.r-project.org")
require(car)
round(vif(Prostate.lm.1),2)
Com podem observar, hi ha diversos valors del FIV que superen el valor 1, per tant reafirmem que tenim un problema de multicol·linealitat.
Per últim calculem el número de condició:
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,"-")
round(max(abs(apply(xc,2,mean))),6)
xs<-apply(xc,2,sd)
x0<-sweep(xc,2,xs,"/")
str(x0)
num [1:97, 1:8] -1.637 -1.989 -1.579 -2.167 -0.508 ... - attr(*, "dimnames")=List of 2 ..$ : chr [1:97] "1" "2" "3" "4" ... ..$ : chr [1:8] "lcavol" "lweight" "age" "lbph" ...
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
y<-Prostate$lpsa
#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(6,-4,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
Prostate.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 coeficients del model ajustat respecte a:
options(repr.plot.width=20,repr.plot.height=7)
old.par<-par(mfrow=c(1,3))
plot(Prostate.ridge.01,xvar="norm")
abline(h=0,lty=2,col='black')
title('Ridge Regression\n')
plot(Prostate.ridge.01,xvar="lambda")
abline(h=0,lty=2,col='black')
title('Ridge Regression\n')
plot(Prostate.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(23025)
Prostate.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<-Prostate.ridge.01.cv$lambda.min
lam1se.ridge1<-Prostate.ridge.01.cv$lambda.1se
options(repr.plot.width=7, repr.plot.height=7)
plot(Prostate.ridge.01.cv)
abline(v=log(bestlam.ridge1),lwd=3,col="cyan")
#Valor de la lambda
round(bestlam.ridge1,2)
#Valor del logaritme de lambda
round(log(bestlam.ridge1),2)
#Valor de 1se
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(Prostate.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.
Tornem a definir un nou ventall de possible valors de $\lambda$:
grid<-exp(seq(0,-4,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.
Prostate.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(Prostate.lasso.01,xvar="norm")
abline(h=0,lty=2,col='black')
title('Lasso Regression\n')
plot(Prostate.lasso.01,xvar="lambda")
abline(h=0,lty=2,col='black')
title('Lasso Regression\n')
plot(Prostate.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(23025)
Prostate.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<-Prostate.lasso.01.cv$lambda.min
lam1se.lasso1<-Prostate.lasso.01.cv$lambda.1se
options(repr.plot.width=7, repr.plot.height=7)
plot(Prostate.lasso.01.cv)
abline(v=log(bestlam.lasso1),lwd=3,col="cyan")
#Valor de la lambda
round(bestlam.lasso1,2)
#Valor del logaritme de lambda
round(log(bestlam.lasso1),2)
#Valor de 1se
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(Prostate.lasso.01,s=bestlam.lasso1,type="coefficients"))
round(coeff.01.lasso,2)
Veiem que 2 dels coeficients de la regressió són exactament iguals a zero, per tant les variables asssociades 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(5,-10,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
Prostate.net.01<-glmnet(x,y,alpha=0.001,lambda=grid)
Prostate.net.02<-glmnet(x,y,alpha=0.3,lambda=grid)
Prostate.net.03<-glmnet(x,y,alpha=0.5,lambda=grid)
Prostate.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(Prostate.net.01,xvar="norm")
abline(h=0,lty=2,col='black')
title('alpha=0.001\n')
plot(Prostate.net.02,xvar="norm")
abline(h=0,lty=2,col='black')
title('alpha=0.3\n')
plot(Prostate.net.03,xvar="norm")
abline(h=0,lty=2,col='black')
title('alpha=0.5\n')
plot(Prostate.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)
Prostate.net.01.cv<-cv.glmnet(x,y,alpha=0.001,lambda=grid)
Prostate.net.02.cv<-cv.glmnet(x,y,alpha=0.3,lambda=grid)
Prostate.net.03.cv<-cv.glmnet(x,y,alpha=0.5,lambda=grid)
Prostate.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<-Prostate.net.01.cv$lambda.min
lam1se.net1<-Prostate.net.01.cv$lambda.1se
plot(Prostate.net.01.cv)
abline(v=log(bestlam.net1),lwd=3,col="cyan")
round(log(bestlam.net1),2)
bestlam.net2<-Prostate.net.02.cv$lambda.min
lam1se.net2<-Prostate.net.02.cv$lambda.1se
plot(Prostate.net.02.cv)
abline(v=log(bestlam.net2),lwd=3,col="cyan")
round(log(bestlam.net2),2)
bestlam.net3<-Prostate.net.03.cv$lambda.min
lam1se.net3<-Prostate.net.03.cv$lambda.1se
plot(Prostate.net.03.cv)
abline(v=log(bestlam.net3),lwd=3,col="cyan")
round(log(bestlam.net3),2)
bestlam.net4<-Prostate.net.04.cv$lambda.min
lam1se.net4<-Prostate.net.04.cv$lambda.1se
plot(Prostate.net.04.cv)
abline(v=log(bestlam.net4),lwd=3,col="cyan")
round(log(bestlam.net4),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.net<-as.numeric(predict(Prostate.net.01,s=bestlam.net1,type="coefficients"))
coeff.02.net<-as.numeric(predict(Prostate.net.02,s=bestlam.net2,type="coefficients"))
coeff.03.net<-as.numeric(predict(Prostate.net.03,s=bestlam.net3,type="coefficients"))
coeff.04.net<-as.numeric(predict(Prostate.net.04,s=bestlam.net4,type="coefficients"))
round(coeff.01.net,2)
round(coeff.02.net,2)
round(coeff.03.net,2)
round(coeff.04.net,2)
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(Prostate.ridge.01)
abline(h=0,lty=2,col='black')
title('alpha=0 (Ridge)\n')
plot(Prostate.net.01,xvar="norm")
abline(h=0,lty=2,col='black')
title('alpha=0.001\n')
plot(Prostate.net.02,xvar="norm")
abline(h=0,lty=2,col='black')
title('alpha=0.3\n')
plot(Prostate.net.03,xvar="norm")
abline(h=0,lty=2,col='black')
title('alpha=0.5\n')
plot(Prostate.net.04,xvar="norm")
abline(h=0,lty=2,col='black')
title('alpha=0.999\n')
plot(Prostate.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(Prostate.ridge.01,s=bestlam.lasso1,type="coefficients"))
round(coeff.01.ridge,5)
print('Coeficients regressió Elastic Net')
coeff.01.net<-as.numeric(predict(Prostate.net.01,s=bestlam.lasso1,type="coefficients"))
coeff.02.net<-as.numeric(predict(Prostate.net.02,s=bestlam.lasso1,type="coefficients"))
coeff.03.net<-as.numeric(predict(Prostate.net.03,s=bestlam.lasso1,type="coefficients"))
coeff.04.net<-as.numeric(predict(Prostate.net.04,s=bestlam.lasso1,type="coefficients"))
round(coeff.01.net,5)
round(coeff.02.net,5)
round(coeff.03.net,5)
round(coeff.04.net,5)
print('Coeficients regressió Lasso')
coeff.01.lasso<-as.numeric(predict(Prostate.lasso.01,s=bestlam.lasso1,type="coefficients"))
round(coeff.01.lasso,5)
[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.