lu2004

Brain Aging Study of Lu et al. (2004)

Descripció

Dades gèniques (403 gens per a 30 mostres) de l’estudi de microarrays de Lu et al. (2004).

Format

  1. 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.

  2. lu2004\$y: edat de cada pacient.

Font

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).

References

Lu, T., et al. (2004), "Gene regulation and DNA damage in the ageing human brain". Nature 429:883–891.

Zuber, V., and K. Strimmer (2011), "High-dimensional regression and variable selection using CAR scores". Statist. Appl. Genet. Mol. Biol. 10 Iss. 1, Article 34.

Objectiu

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.

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.

Mínims quadrats ordinaris (MQO)

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.

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.

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.

Ridge

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.

Primer definim un ventall de possible valors de $\lambda$:

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

Ara utilitzem la funció plotper realitzar el gràfic que ens permet veure l'estimació dels coeficientes del model ajustat respecte a:

  1. La norma L1
  2. El log($\lambda$)
  3. El percentatge de desviació explicada

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.

Tornem a utilitzar la funció plot per veure gràficament els resultats de la validació encreuada. Podem veure dues marques:

  1. La blau marca el log($\lambda_{min}$) on $\lambda_{min}$ proporciona un error mitjà de la validació encreuada mínim.
  2. La negra marca 1se és el valor més gran de $\lambda$ de manera que es troba a 1 desviament estàndard de la $\lambda_{min}$

Finalment utilitzem la funció predict() per reajustar els coeficients utilitzant el valor de la $\lambda$ resultant de la validació encreuada.

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.

Lasso

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$:

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.

Ara utilitzem la funció plotper realitzar el gràfic que ens permet veure l'estimació dels coeficientes del model ajustat respecte a:

  1. La norma L1
  2. El log($\lambda$)
  3. El percentatge de desviació explicada

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.

Tornem a utilitzar la funció plot per veure gràficament els resultats de la validació encreuada. Podem veure dues marques:

  1. La blau marca el log($\lambda_{min}$) on $\lambda_{min}$ proporciona un error mitjà de la validació encreuada mínim.
  2. La negra marca 1se és el valor més gran de $\lambda$ de manera que es troba a 1 desviament estàndard de la $\lambda_{min}$

Finalment utilitzem la funció predict() per reajustar els coeficients utilitzant el valor de la $\lambda$ resultant de la validació encreuada.

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.

Elastic Net

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$:

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

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:

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.

Tornem a utilitzar la funció plot per veure gràficament els resultats de la validació encreuada. Podem veure dues marques:

  1. La blau marca el log($\lambda_{min}$) on $\lambda_{min}$ proporciona un error mitjà de la validació encreuada mínim.
  2. La negra marca 1se és el valor més gran de $\lambda$ de manera que es troba a 1 desviament estàndard de la $\lambda_{min}$

Finalment utilitzem la funció predict() per reajustar els coeficients utilitzant el valor de la $\lambda$ resultant de la validació encreuada.

Similars als coeficients de la regressió lasso, amb alguns coeficients nuls, que resolen les limitacions del lasso.

Comparacions

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.

Finalment, podem comparar els coeficients obtinguts mitjançant la lambda que minimitza la mitjana de l'error per a la 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.