Técnicas de selección de variables en minería estadística de datos

FACULTAD DE MATEM´ ATICAS ESTAD´ ISTICA E INVESTIGACI´ ON OPERATIVA Trabajo Fin de Grado T´ecnicas de Selecci´on de Variables en Miner´ıa Estad´ıstica de Datos Adri´an Guerra de la Corte Dirigido por: D˜na. Inmaculada Barranco Chamorro Sevilla, Junio 2016. Abstract A common problem in data mining, when statistical regression models are used, is to choose properly the variables to be included in the model. Throughout this work the main statistical techniques for the selection and regularization of variables will be reviewed. Also applications of these techniques will be performed by using R. The work is divided into four chapters. In Chapter 1, we review the linear regression model, and the different correlation coefficients. In this way we introduce the basic tools to study methods of selection and regularization of variables in linear regression models. In Chapter 2, we will see the most common criteria used for the selection of variables in classical linear models. So, we will deal with: Adjusted coefficient of determination,Mallow’s Coefficient,Cross Validation method,Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC). These criteria will be compared between them. Also, the main problems we may have in practice when using multiple linear regression techniques are studied. An application in Rhas been included to illustrate the performance of the different methods. In Chapter 3, we focus on the so-called heuristic methods, which are a first approach to the problem of selection of variables when we have a very large number of regressors. So, selection techniques such as forward, backward and step by step are studied. Their use is again illustrated with an application. In Chapter 4, we discuss the regularization techniques. We focus on ridge regression and LASSO regression. In this context, we show that by applying regularization techniques the problem becomes manageable, since a set of restrictions is imposed on the set of admissible solutions. As well, the geometric properties of the estimators are studied. As before, an application is included to illustrate the use of the discussed techniques in the field of medecine. Finally, the work is completed by an appendix, which contains the Rand Mathematica codes implemented for the development of the figures, as well as the packages of Rused, and the literature consulted. III IV Resumen Al utilizar modelos de regresi´on en Miner´ıa Estad´ıstica de Datos, un problema com´un es elegir de forma adecuada las variables a incluir en el modelo. A lo largo de este trabajo se revisar´an las t´ecnicas estad´ısticas que existen para la selecci´on y regularizaci´on de variables. As´ı mismo se realizar´an aplicaciones de dichas t´ecnicas, b´asicamente con el software R. El trabajo se estructura en cuatro cap´ıtulos. En el Cap´ıtulo 1, revisamos el modelo de regresi´on lineal, as´ı como los diferentes coeficientes de correlaci´on. De esta forma introducimos las herramientas b´asicas para abordar el estudio de los m´etodos de selecci´on y regularizaci´on de variables en los modelos de regresi´on lineal. En el Cap´ıtulo 2, veremos los criterios m´as usados para la selecci´on de variables en modelos lineales cl´asicos. Se recogen as´ı: el coeficiente de determinaci´on corregido o ajustado, el coeficiente Cpde Mallows, el m´etodo de validaci´on cruzada, el criterio de informaci´on de Akaike (AIC) y el criterio de informaci´on bayesiana (BIC). Se realizan comparaciones entre ellos, y se recogen los principales problemas que se nos pueden presentar en la pr´actica al utilizar las t´ecnicas de regresi´on lineal m´ultiple. As´ı mismo, cabe destacar que se ha ilustrado el uso de las distintas t´ecnicas expuestas con una aplicaci´on realizada con R. En el Cap´ıtulo 3, nos centraremos en los llamados m´etodos heur´ısticos, los cuales son una primera aproximaci´on al problema de selecci´on de variables cuando tenemos un n´umero muy grande de variables regresoras. Se recogen las denominadas t´ecnicas de selecci´on hacia adelante, hacia atr´as y paso a paso. Su uso se ilustra de nuevo con una aplicaci´on. En el Cap´ıtulo 4, trataremos las t´ecnicas de regularizaci´on, principalmente el modelo de regresi´on contra´ıda (ridge regression) y el modelo de regresi´on LASSO (LASSO regression). Estas t´ecnicas permiten solventar las dificultades que surgen cuando se presentan problemas de colinealidad o soluciones num´ericas inestables. En este contexto, mostramos que regularizar significa, hacer el problema tratable, imponiendo una serie de restricciones al conjunto de soluciones admisibles. Adem´as se estudian las propiedades geom´etricas de V los estimadores obtenidos. De nuevo se incluye una aplicaci´on, en el campo de la Medicina, que ilustra el uso de las t´ecnicas expuestas. Finalmente, el trabajo se completa con un anexo, en el que se recogen los c´odigos Ry de Mathematica implementados para la elaboraci´on de las figuras, as´ı como los paquetes de Rutilizados, y la bibliograf´ıa consultada. VI ´ Indice general Abstract III Resumen V 1. Conceptos previos 1 1.1. Modelo de regresi´on lineal simple . . . . . . . . . . . . . . . . 1 1.1.1. C´alculo de los estimadores . . . . . . . . . . . . . . . . 2 1.2. Modelo de regresi´on lineal m´ultiple . . . . . . . . . . . . . . . 4 1.2.1. C´alculo de los estimadores . . . . . . . . . . . . . . . . 5 1.2.2. Intepretaci´on de los coeficientes en un modelo de regresi´on lineal m´ultiple . . . . . . . . . . . . . . . . . . 8 1.2.3. Contrastes......................... 8 1.3. Coeficientes de correlaci´on . . . . . . . . . . . . . . . . . . . . 9 1.3.1. Relaci´on entre las correlaciones parciales y la m´ultiple . 12 2. T´ecnicas de selecci´on de variables en modelos lineales cl´asicos 15 2.1. Coeficiente de determinaci´on corregido o ajustado . . . . . . . 16 2.1.1. Aplicaci´on......................... 17 2.2. Coeficiente CpdeMallows .................... 22 2.2.1. Aplicaci´on......................... 23 2.3. Validaci´on cruzada . . . . . . . . . . . . . . . . . . . . . . . . 25 2.3.1. Validaci´on cruzada en r iteraciones . . . . . . . . . . . 25 2.3.2. Validaci´on crudada dejando uno fuera . . . . . . . . . . 25 2.3.3. Aplicaci´on......................... 26 2.4. Criterio de Informaci´on de Akaike . . . . . . . . . . . . . . . . 29 2.4.1. Aplicaci´on......................... 31 2.5. Criterio de Informaci´on Bayesiana . . . . . . . . . . . . . . . . 33 2.5.1. Aplicaci´on......................... 34 2.6. Comparaci´on de criterios . . . . . . . . . . . . . . . . . . . . . 35 2.7. Problemas en la regresi´on m´ultiple . . . . . . . . . . . . . . . 36 VII 2.7.1. Error de especificaci´on . . . . . . . . . . . . . . . . . . 36 2.7.2. Hip´otesis de normalidad . . . . . . . . . . . . . . . . . 37 2.7.3. Robustez.......................... 39 2.7.4. Heterocedasticidad . . . . . . . . . . . . . . . . . . . . 41 2.7.5. Multicolinealidad . . . . . . . . . . . . . . . . . . . . . 44 3. M´etodos heur´ısticos para la selecci´on de variables 51 3.1. Selecci´on hacia delante . . . . . . . . . . . . . . . . . . . . . . 52 3.1.1. Aplicaci´on......................... 52 3.2. Selecci´on hacia atr´as . . . . . . . . . . . . . . . . . . . . . . . 54 3.2.1. Aplicaci´on......................... 55 3.3. Selecci´on paso a paso . . . . . . . . . . . . . . . . . . . . . . . 56 3.3.1. Aplicaci´on......................... 56 4. T´ecnicas de regularizaci´on 59 4.1. Regresi´on contra´ıda . . . . . . . . . . . . . . . . . . . . . . . . 60 4.1.1. Aplicaci´on......................... 62 4.2. Regresi´on LASSO . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.2.1. Aplicaci´on......................... 68 4.3. Propiedades geom´etricas de los estimadores regularizados . . . 71 A. Anexo 76 A.1. Comandos en Rdelasgr´aficas.................. 76 A.1.1.Figura2.1......................... 76 A.1.2.Figura2.2......................... 76 A.1.3.Figura4.1......................... 76 A.2. Comandos en Mathematica de las gr´aficas . . . . . . . . . . . 77 A.2.1.Figura4.3......................... 77 A.3. Paquetes de R........................... 78 VIII Cap´ıtulo 1 Conceptos previos En este cap´ıtulo explicaremos los resultados b´asicos a la hora de introducir y comprender el estudio de m´etodos para la selecci´on adecuada de las variables a incluir en un modelo de regresi´on lineal, simple y m´ultiple, as´ı como las t´ecnicas de regularizaci´on de dichas variables. Dividiremos el cap´ıtulo en 3 secciones: regresi´on lineal simple, regresi´on lineal m´ultiple y coeficientes de correlaci´on. Esta ´ultima secci´on tendr´a especial relevancia en el siguiente cap´ıtulo. 1.1. Modelo de regresi´on lineal simple Comenzaremos con el modelo de regresi´on lineal simple, que consiste en expresar la dependencia lineal de la variable objetivo o dependiente, y, respecto a otras dos variables: la variable independiente, explicativa o covariable, x, y el t´ermino error o perturbaci´on del modelo, uas´ı yi=β0+β1xi+ui,con (xi, yi) variables num´ericas donde yiyuison variables aleatorias, xies una variable conocida una vez observada yi, y β0yβ1son par´ametros desconocidos del modelo. Las hip´otesis del modelo pueden formularse en t´erminos de la variable perturbaci´on, ui, o de forma equivalente en t´erminos de la variable dependiente, y. As´ı podemos establecer las siguientes hip´otesis: La perturbaci´on debe tener esperanza nula, es decir E(ui) = 0 ⇔E(yi) = β0+β1xi. 1 1.2.2. Intepretaci´on de los coeficientes en un modelo de regresi´on lineal m´ultiple El coeficiente de regresi´on estimado para una variable xien el modelo de regresi´on lineal m´ultiple, ˆ βi, representa el efecto sobre la variable objetivo cuando la variable xiaumenta en una unidad y las dem´as variables explicativas premanecen constantes. Puede interpretarse como el efecto diferencial de esta variable cuando eliminamos o controlamos los efectos de las dem´as. Debemos distinguir dos situaciones a la hora de interpretar los coeficientes, cuando las variables explicativas est´an incorreladas y cuando no lo est´an: Cuando todas las variables explicativas est´an incorreladas se calcula de la misma manera que en la regresi´on simple. Pues en este caso, el efecto diferencial de la variable, medido por la regresi´on m´ultiple, es igual al efecto total medido por la regresi´on simple. Cuando las variables est´an correladas, el coeficiente de regresi´on de xi se puede expresar tambi´en como el cociente entre una covarianza y una varianza. Con la salvedad de que en la regresi´on simple se utiliza la covarianza entre la variable objetivo y xi, y en la m´ultiple se utiliza la covarianza entre la variable objetivo y la parte diferencial de xio no correlada con el resto de variables explicativas. La parte diferencial de xiest´a definida por los residuos de una regresi´on entre la variable xi y el resto de variables explicativas en la ecuaci´on de regresi´on. Estos residuos ser´an ei,R. Luego, ˆ βi=Cov(y, ei,R)/V ar(ei,R), donde no se usa la variable xicomo en la regresi´on simple, sino la parte diferencial de ella, ei,R. En el caso de que esta variable s´ı est´e incorrelada, xi=ei,R y el coeficiente de la regresi´on m´ultiple es igual al de la regresi´on simple. 1.2.3. Contrastes A la hora de calcular un modelo de regresi´on lineal, los contrastes son una herramienta importante. En esta secci´on hablaremos de los dos contrastes que usaremos en este trabajo, el contraste global de regresi´on y el contraste individual de la t. Contraste global de regresi´on El contraste es el siguiente: H0:β1=... =βk= 0. 8 H1: alg´un βi6= 0, i = 1, ..., k. Y el estad´ıstico resultante se denota por Fy se calcula F=SCE/k SCR/n −k−1=ˆs2 expl ˆs2 r . Bajo la hip´otesis H0,F∼Fk,n−k−1. Dicho contraste se traduce en Si acepto H0⇒ninguna de las variables explicativas consideradas influyen linealmente en la variable respuesta. Si rechazo H0⇒alguna o todas las variables explicativas consideradas influyen linealmente en la variable respuesta. Contraste individual de la t Para cada variable se plantea el siguiente contraste: H0:βi= 0. H1:βi6= 0. El estad´ıstico que resulta del contraste est´a basado en el estad´ıstico de Wald y se define como sigue ti=ˆ βi bsr√qii siendo ˆ βiel estad´ıstico de Wald yqii el t´ermino (ii) de la matriz (X0X)−1. ´ Este, bajo H0, sigue una distribuci´on tn−k−1. Como es usual, en contrastes de hip´otesis, se rechaza H0, si el p-valor obtenido es menor o igual que el nivel de significaci´on del contraste, α. En este caso podremos suponer que βi6= 0. Si no se rechaza H0, podremos suponer que βi= 0. 1.3. Coeficientes de correlaci´on En esta secci´on se definen y estudian las propiedades de los coeficientes de correlaci´on lineal simple, coeficiente de determinaci´on, coeficiente de correlaci´on m´ultiple yparcial, as´ı como las relaciones existentes entre ellos. Una medida de la relaci´on lineal entre dos variables cualesquiera es el coeficiente de correlaci´on lineal simple. 9 Definici´on 1.3.1. Dadas dos variables xey, se denomina coeficiente de correlaci´on lineal simple,rxy, a rxy =Cov(x, y) sxsy donde sxysyson las desviaciones t´ıpicas muestrales de las variables xey, respectivamente. Dicho coeficiente se puede expresar en funci´on de la varianza residual r2 xy =SCE(x, y) SCT(y)= 1 −SCR(x, y) SCT(y). Definici´on 1.3.2. Definimos el coeficiente de determinaci´on,R2, de un modelo, para evaluar la bondad de ajuste de una recta de regresi´on (simple o m´utiple) con una proporci´on de la variaci´on explicada con la siguiente expresi´on R2=SCE SCT =P(ˆyi−y)2 P(yi−y)2,0≤R2≤1. Siendo √R2el coeficiente de correlaci´on m´ultiple. Definici´on 1.3.3. Dado un conjunto de variables explicativas (x1, ..., xk), el coeficiente de correlaci´on parcial, denotado por rij,1,...,i−1,i+1,...,j−1,j+1,...,k, entre dos cualesquiera de ellas, xiyxj, mide la relaci´on lineal entre xiyxj una vez eliminados los efectos de las dem´as sobre ellas. Se denotar´a por r12,34...k. La manera de calcularlo es ´utilizando la siguiente expresi´on: e1,34..k =ˆ βe2,34..k +u donde e1,34..k ye2,34..k son los residuos de la regresi´on m´ultiple de x1yx2respecto a las dem´as variables (x3, ..., xk), de esta manera obtendr´ıamos r12,34..k. C´alculo del coeficiente de correlaci´on parcial A continuaci´on se deduce la f´ormula de la correlaci´on parcial de dos variables (x, y) cuando se mantiene constante una tercera variable z. 10 Proposici´on 1.3.1. El coeficiente de correlaci´on parcial entre xyymanteniendo constante z, viene dado por rxy.z =rxy −rxzryz √(1−r2 xz)(1−r2 yz).(1.3) donde rxy, ..., ryz son los coeficientes de correlaci´on lineal simple entre las variables implicadas. Demostraci´on. Supondremos que las tres variables tienen media cero para simplicar la exposici´on, esto no altera el resultado. Sean las rectas ˆx=az, ˆy=bz cuyos coeficientes son, por definici´on: a=Pxizi Pz2 i b=Pyizi Pz2 i . Puesto que la correlaci´on parcial entre xey, fijada z, es la correlaci´on entre los residuos de estas regresiones, tenemos que: Correlaci´on[(x−ˆx)(y−ˆy)] = Cov(x, y) pV ar(x−ˆx)V ar(y−ˆy)=rxy.z. Pasemos a calcular el numerador: n Cov(x, y) = X(xi−ˆxi)(yi−ˆyi) = X(xi−azi)(yi−bzi) =Xxiyi−aXziyi−bXzixi+ab Xz2 i. Sustituyendo ayben la expresi´on anterior, llegamos a: n Cov(x, y) = Xxiyi−(Pxizi)(Pziyi) Pz2 i−(Pyizi)(Pxizi) Pz2 i +(Pxizi)(Pyizi) Pz2 i =Xxiyi−(Pxizi)(Pziyi) Pz2 i . Si introducimos ahora los coeficientes de correlaci´on simples, nos queda: n Cov[(x−ˆx)(y−ˆy)] = rxyqXx2 iXy2 i−rxzryzqXx2 iXy2 i = (rxy −rxzryz)qXx2 iXy2 i. 11 Con eso tenemos calculado el numerador, pasemos a hallar denominador: n V ar(x−ˆx) = X(xi−azi)2=Xx2 i−(Pxizi)2 Pz2 i n V ar(y−ˆy) = X(yi−bzi)2=Xy2 i−(Pyizi)2 Pz2 i . Sustituyendo en las varianzas de los residuos los coeficientes de correlaci´on simples, se tiene: n V ar(x−ˆx) = Xx2 i(1 −r2 xz) n V ar(y−ˆy) = Xy2 i(1 −r2 yz). Por lo que finalmente llegamos a la expresi´on 1.3.  Este coeficiente al cuadrado tiene la misma interpretaci´on que el coeficiente de correlaci´on simple. Es decir, r2 xy.z representa la proporci´on de variaci´on explicada respecto a la variaci´on no explicada por otra regresi´on previa. 1.3.1. Relaci´on entre las correlaciones parciales y la m´ultiple Supongamos, por simplificar, que tenemos nada m´as dos variables explicativas, x1yx2. Sea ryx1el coeficiente de correlaci´on simple entre la variable objetivo y x1. Que por lo visto anteriormente es: r2 yx1=SCE(y, x1) SCT(y)= 1 −SCR(y, x1) SCT(y) donde SCE(y, x1) es la variaci´on explicada en la regresi´on de yrespecto a x1. Luego: SCR(y, x1) = SCT(y)(1 −r2 yx1). Ahora debemos determinar la parte diferencial de la segunda variable, e2,1, los cuales calculamos haciendo la regresi´on x2=ˆ bx1+e2,1. Una vez calculados, los relacionamos con la parte de la variable objetivo que no est´a explicada por x1, que ser´an los residuos ey.x1de la regresi´on simple de yrespecto a x1. La relaci´on de ambos residuos est´a dada por el coeficiente de correlaci´on parcial, el cual nos proporciona los residuos de la regresi´on m´ultiple. La estimaci´on ser´a ey.x1=ˆ βe2,1+ey,12, por tanto: r2 y2,1= 1 −SCR(ey,12) SCT(ey.x1)= 1 −SCR(ey,12) SCR(y, x1). 12 Y como ey,12 son los residuos de la regresi´on m´ultiple con ambas variables: R2= 1 −SCR(ey,12) SCT(y) luego 1−R2= (1 −r2 yx1)(1 −r2 y2,1). Que se interpreta como la proporci´on de la variabilidad no explicada en la regresi´on m´ultiple es el producto de: la proporci´on no explicada en la regresi´on simple de la variable objetivo yx1. la proporci´on no explicada en la regresi´on de la variable objetivo y x2 con x1fija. Dicho resultado se puede extender para kregresores: 1−R2= (1 −r2 y1)(1 −r2 y2,1)(1 −r2 y3,12)...(1 −r2 yk,12...k−1). Con esta idea tambi´en podemos relacionar los coeficientes de correlaci´on m´ultiple con kyk−1 variables, en funci´on del coeficiente de correlaci´on parcial de la variable no incluida, xh. Llamando R2 kyR2 k−1a los coeficientes de correlaci´on m´ultiple con kyk−1 variables, respectivamente, llegamos a 1−R2 k= (1 −R2 k−1)(1 −r2 yh,12..k). lo que es lo mismo que: R2 k−R2 k−1=r2 yh,12..k(1 −R2 k−1). El t´ermino de la izquierda representa el incremento de variaci´on explicada entre la regresi´on que incluye a la variable y la que no la incluye. El t´ermino de la derecha es el producto del porcentaje de variaci´on explicada por xhrespecto a la variaci´on no explicada por las restantes variables y del porcentaje de variaci´on no explicada por las restantes variables, x1, ..., xkrespecto al total. Esta expresi´on nos permite calcular los coeficientes de correlaci´on parcial de cada variable a partir de un programa que nos calcule la regresi´on. En general, la f´ormula anterior puede escribirse de tal manera que, notando (¯ h), por el modelo que no incluye a xh, y por (h) al modelo que s´ı la incluye: ∆SCE(h) SCT =∆SCE(h) SCR(¯ h)·SCR(¯ h) SCT 13 donde ∆SCE(h) = SCE(todas)−SCE(¯ h). Por ´ultimo, utilizando el estad´ıstico tmencionado en la Secci´on 1.2.3 para contrastar la hip´otesis de βh= 0, llegamos a: r2 yh,12...k =t2 h t2 h+n−(k+ 1) lo cual nos permite calcular el coeficiente de correlaci´on parcial si sabemos el estad´ıstico tpara ese coeficiente. 14 Cap´ıtulo 2 T´ecnicas de selecci´on de variables en modelos lineales cl´asicos A la hora de construir un modelo tenemos diferentes posibilidades, las cuales se ajustan mejor o peor a la realidad. En este cap´ıtulo nos centraremos en los criterios m´as usados para la selecci´on de variables en modelos lineales, que son: Coeficiente de determinaci´on corregido o ajustado: Es un coeficiente que mide la intensidad de la relaci´on lineal entre la variable objetivo y las predictoras. Coeficiente Cpde Mallows: Criterio que recibe el nombre del estad´ıstico brit´anico Colin Lingwood Mallows. Este criterio selecciona el modelo que tiene mayor capacidad de predicci´on en vez del que est´a mejor ajustado. La capacidad de predicci´on se mide con el error cuadr´atico medio (ECM). Validaci´on cruzada: Evoluci´on del llamado holdout method que se basa en la partici´on del conjunto de datos en dos, uno nos permite estimar los par´ametros del modelo, y el otro evaluar la capacidad predictiva de ´este. De esta forma se selecciona el modelo valorando su bondad de ajuste y capacidad de predicci´on. Criterio de Informacion de Akaike (AIC): Criterio propuesto por el estad´ıstico japon´es Hirotugu Akaike y que est´a basado en la teor´ıa de la informaci´on. Est´a definido de forma que bonifica la bondad de ajuste y penaliza la inclusi´on de par´ametros a estimar, lo que ayuda a evitar el fen´omeno del sobreajuste. 15 Criterio de Informaci´on Bayesiana (BIC): El profesor Gideon E.Schwarz propuso este criterio bajo un enfoque bayesiano que se basa en las probabilidades a posteriori de los modelos. Es, junto al AIC, el m´as usado. 2.1. Coeficiente de determinaci´on corregido o ajustado El coeficiente de determinaci´on R2es una medida de bondad de ajuste de un modelo a unos datos. Recordemos que R2nos da la proporci´on de la variabilidad de Yaplicada por el modelo, es decir, R2=SCExplicada SCT =Pn i=1(ˆyi−y)2 Pn i=1(yi−y)2= 1 −SCResidual SCT y adem´as 0 ⩽R2⩽1. Cuanto m´as cercano est´e a 1, mejor ajustado est´a el modelo. R2no nos sirve para comparar modelos diferentes, puesto que siempre aumenta cuando se a˜naden nuevas variables explicativas al modelo, lo que nos llevar´ıa a tomar modelos con innumerables variables superfluas. Por eso surge el coeficiente de determinaci´on corregido o ajustado, que sirve para solventar este problema puesto que incluye un t´ermino de correcci´on por el n´umero de par´ametros en el modelo. b R2 aj,k = 1 −n−1 n−k(1 −R2 k). b R2 aj es muy popular y viene incorporado en los programas estad´ısticos y resulta de especial inter´es en situaciones en las que el n´umero de variables explicativas est´a cercano al n´umero de observaciones de la muestra. Teorema 2.1.1. El estad´ıstico b R2 aj,k aumenta al introducir un nuevo par´ametro, βk+1, en la ecuaci´on de regresi´on si el estad´ıstico Qhasociado al contraste de significaci´on de dicho par´ametro es mayor que 1. Qhse define como: Qh=SCRk−SCRk+1 SCRk+1 ×n−k−1 1 donde SCRkes la suma de cuadrados de los residuos en el modelo con k covariables. 16 Demostraci´on. Para hacer el contraste del (k+1)-´esimo par´ametro emplearemos el estad´ıstico Qh, que se defini´o como: Qh=SCRk−SCRk+1 SCRk+1 ×n−k−1 1 =R2 k+1 −R2 k 1−R2 k+1 ×n−k−1 1 por tanto: (1 −R2 k+1)Qh= (R2 k+1 −R2 p)(n−k−1) Qh−QhR2 k+1 = (n−k−1)R2 k+1 −(n−k−1)R2 k Qh+ (n−k−1)R2 k=R2 k+1[(n−k−1) + Qh] despejando R2 k+1: R2 k+1 =Qh+ (n−k−1)R2 k (n−k−1) + Qh = 1 n−k−1Qh+R2 k 1 + 1 n−k−1Qh . Sustituyendo esta expresi´on en la definici´on de b R2 aj,k+1, tenemos: b R2 aj,k+1 = 1 −(1 −R2 k+1)n−1 n−k−1= 1 −1−R2 k n−k−1+Qh n−k−1×n−1 n−k−1 = 1 −(1 −R2 k)n−1 n−k−1 + Qh = 1 −(1 −R2 k)n−1 n−k | {z } b R2 aj,k ×n−k n−k−1 + Qh | {z } t de lo que se deduce que b R2 aj,k+1 ≥b R2 ksi Qh>1.  2.1.1. Aplicaci´on A continuaci´on se recoge un ejemplo que ilustra el uso del coeficiente de determinaci´on corregido o ajustado. Utilizaremos el conjunto de datos que nos proporciona Fahrmeir, L. et al. [1]. En primer lugar se cargar´a el conjunto de datos. > golf <- read.table("golffull.txt", header=TRUE) > attach(golf) 17 Call: mle.cp(formula = price ~ kilometer + age + extras1 + extras2 + TIA, data = golf) Mallows Cp: (Intercept) kilometer age extras1 extras2 TIA cp [1,] 1 1 1 1 0 0 2.422 [2,] 1 1 1 1 0 1 4.005 [3,] 1 1 1 1 1 0 4.407 [4,] 1 1 1 1 1 1 6.000 Printed the first 4 best models En nuestro caso, como tenemos que seleccionar los modelos cuyas variables age ykilometer est´en presentes, hemos tenido que hacer alguna modificaci´on > misres<- subset(cp$cp, cp$cp[,2]==1 & cp$cp[,3]==1) > head(misres) (Intercept) kilometer age extras1 extras2 TIA cp [1,] 0 1 1 0 0 0 622.527918 [2,] 1 1 1 0 0 0 3.551917 [3,] 0 1 1 1 0 0 590.044015 [4,] 1 1 1 1 0 0 2.422241 [5,] 0 1 1 0 1 0 608.241325 [6,] 1 1 1 0 1 0 5.505495 Vemos que los modelos que manejamos contienen a las variables deseadas. Ahora debemos escoger el que tenga menor cp, con lo que vamos a ordenarlos: > ordenado<-misres[ order(misres[,7]), ] > head(ordenado) (Intercept) kilometer age extras1 extras2 TIA cp [1,] 1 1 1 1 0 0 2.422241 [2,] 1 1 1 0 0 0 3.551917 [3,] 1 1 1 1 0 1 4.005040 [4,] 1 1 1 1 1 0 4.406582 [5,] 1 1 1 0 0 1 5.401711 [6,] 1 1 1 0 1 0 5.505495 En esta tabla se muestran los 6 mejores modelos, seg´un este criterio, que tienen las variables explicativas impuestas por nosotros. Luego el mejor modelo ser´a: price ∼intercept +kilometer +age +extras1. 24 2.3. Validaci´on cruzada El primer paso y com´un a todos los tipos de validaciones cruzadas es dividir el conjunto de datos que tenemos en dos tipos de conjuntos: Un tipo de conjunto que nos sirve para estimar los par´ametros del modelo, llamado training set. Un tipo de conjunto de validaci´on que sirve para valorar la capacidad predictiva del modelo, llamado testing set. Pasemos a ver las dos validaciones cruzadas m´as usadas. 2.3.1. Validaci´on cruzada en r iteraciones La validacion cruzada en r iteraciones (r-fold cross validation) comienza agrupando los datos en rsubconjuntos de tama˜no similar. Uno se utiliza para validar, y los restantes (r−1) se consideran como conjuntos de estimaci´on. Repetiremos este proceso rveces, una vez con cada uno de los subconjuntos de validaci´on. Nos quedamos con el modelo en el que la suma de los errores de predicci´on al cuadrado sea m´as peque˜no, es decir: m´ın{CV }dondeCV =1 n n X i=1 (yi−ˆyiM )2}. Normalmente se suele utilizar el 10-fold cross validation ´o 5-fold cross validation, dependiendo del tama˜no de datos que tengamos. 2.3.2. Validaci´on crudada dejando uno fuera Un caso importante del anterior, que cabe destacar, es la validaci´on cruzada dejando una sola observaci´on fuera (leave-one-out cross validation). En este caso, el error es muy peque˜no, en cambio, el coste computacional es elevado, pues hay que calcular niteraciones y analizar para cada iteraci´on los datos de ambos conjuntos. El estad´ıstico para decidir en este caso es: CV =1 n n X i=1 (yi−ˆy−i iM )2 donde, ˆy−i iM ≡estimaci´on cuando se ha eliminado la observaci´on i-´esima. Se tiene una expresi´on sencilla para este coeficiente, sin tener que rehacer todos los c´alculos, bas´andonos en los ˆyiM originales: CV =1 n n X i=1 yi−ˆyiM 1−hiiM 2 25 donde 1 −hii,M son los elementos diagonales de la matriz hat. Es importante destacar que ambos tipos de validaci´on tienen ciertas limitaciones: El training set y el testing set deben ser extra´ıdos de la misma poblaci´on, en caso de no serlo, la validaci´on no producir´ıa resultados significativos. Esta herramienta no es v´alida cuando tenemos un sistema que evoluciona con el tiempo, pues podr´ıa darse el caso de que ambos conjuntos mencionados anteriormente sufrieran cambios sistem´aticos, por ejemplo: si tenemos un modelo que utilizamos para predecir el valor de las acciones, el cual ha sido calculado en un training set en un periodo de tiempo determinado, ´este no ser´a eficiente a la hora de predecir el valor de la misma poblaci´on en el siguiente periodo de tiempo. Es importante notar que en los datos del training set debemos evitar que haya alg´un dato que est´e tambi´en en el testing set. 2.3.3. Aplicaci´on Para ilustrar este m´etodo procederemos inicialmente de manera similar a como se hizo en 2.1.1 y teniendo en cuenta el estudio previo realizado en ´el, vamos a calcular los distintos modelos para ver c´ual tiene menor CV . Para calcular los distintos CV usaremos el paquete lattice, necesario para el paquete DAAG, el cual contiene a la funci´on CV lm. Dicha funci´on realiza las m-fold cross validation seg´un el valor que le asignemos a m. Comencemos leyendo los datos, definiendo los modelos y cargando los paquetes necesarios: > golf <- read.table("golffull.txt", header=TRUE) > attach(golf) > mod1 <- lm(price~kilometer+age, data=golf) > mod2 <- lm(price~kilometer+age+extras1, data=golf) > mod3 <- lm(price~kilometer+age+extras2, data=golf) > mod4 <- lm(price~kilometer+age+TIA, data=golf) > mod5 <- lm(price~kilometer+age+extras1+extras2, data=golf) > mod6 <- lm(price~kilometer+age+extras1+TIA, data=golf) > mod7 <- lm(price~kilometer+age+extras2+TIA, data=golf) > mod8 <- lm(price~kilometer+age+extras1+extras2+TIA, data=golf) > library(lattice) > library(DAAG) 26 Una vez que hemos escrito los argumentos en la funci´on, ´esta nos muestra la tabla con el an´alisis de la varianza, los resultados de cada fold, y al final la suma de los errores al cuadrado (MS). > RES1<-CVlm(data=golf, form.lm=mod1, m=10) Analysis of Variance Table Response: price Df Sum Sq Mean Sq F value Pr(>F) kilometer 1 88.1 88.1 146 <2e-16 *** age 1 75.2 75.2 124 <2e-16 *** Residuals 169 102.2 0.6 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 fold 1 Observations in test set: 17 9 11 16 45 50 55 62 64 73 Predicted 4.93 4.502 4.557 4.17 2.96 4.12 5.12 4.88 3.756 cvpred 4.92 4.488 4.544 4.16 2.95 4.11 5.12 4.88 3.746 price 6.35 3.823 4.950 2.90 4.20 3.99 6.15 4.50 3.000 CV residual 1.43 -0.665 0.406 -1.26 1.25 -0.12 1.03 -0.38 -0.746 85 101 118 125 148 150 155 164 Predicted 3.5638 4.31 2.827 3.8879 2.210 2.749 3.22 2.721 cvpred 3.5542 4.31 2.817 3.8861 2.202 2.744 3.22 2.721 price 3.6500 3.20 3.800 3.9500 1.900 3.100 4.20 2.400 CV residual 0.0958 -1.11 0.983 0.0639 -0.302 0.356 0.98 -0.321 Sum of squares = 11.1 Mean square = 0.65 n = 17 . . . fold 10 Observations in test set: 17 5 10 13 20 31 32 34 40 43 Predicted 5.19 5.383 4.94 3.88 3.405 3.89 4.445 4.558 5.49 cvpred 5.11 5.274 4.87 3.90 3.444 3.88 4.391 4.489 5.34 27 price 6.20 5.900 5.55 2.50 3.250 2.40 4.600 5.450 7.00 CV residual 1.09 0.626 0.69 -1.40 -0.194 -1.48 0.209 0.961 1.66 44 76 83 110 117 119 133 137 Predicted 3.33 3.690 3.185 3.766 2.7887 2.334 3.0203 2.594 cvpred 3.37 3.676 3.211 3.722 2.8252 2.411 3.0215 2.632 price 3.70 2.800 2.600 3.900 2.9000 1.450 2.9990 1.950 CV residual 0.33 -0.876 -0.611 0.178 0.0748 -0.961 -0.0225 -0.682 Sum of squares = 12.7 Mean square = 0.75 n = 17 Overall (Sum over all 17 folds) ms 0.617 23456 2 3 4 5 6 7 Predicted (fit to all data) price Small symbols show cross−validation predicted values Fold 1 Fold 2 Fold 3 Fold 4 Fold 5 Fold 6 Fold 7 Fold 8 Fold 9 Fold 10 Fold 1 Fold 2 Fold 3 Fold 4 Fold 5 Fold 6 Fold 7 Fold 8 Fold 9 Fold 10 Fold 1 Fold 2 Fold 3 Fold 4 Fold 5 Fold 6 Fold 7 Fold 8 Fold 9 Fold 10 Fold 1 Fold 2 Fold 3 Fold 4 Fold 5 Fold 6 Fold 7 Fold 8 Fold 9 Fold 10 Fold 1 Fold 2 Fold 3 Fold 4 Fold 5 Fold 6 Fold 7 Fold 8 Fold 9 Fold 10 Fold 1 Fold 2 Fold 3 Fold 4 Fold 5 Fold 6 Fold 7 Fold 8 Fold 9 Fold 10 Fold 1 Fold 2 Fold 3 Fold 4 Fold 5 Fold 6 Fold 7 Fold 8 Fold 9 Fold 10 Fold 1 Fold 2 Fold 3 Fold 4 Fold 5 Fold 6 Fold 7 Fold 8 Fold 9 Fold 10 Fold 1 Fold 2 Fold 3 Fold 4 Fold 5 Fold 6 Fold 7 Fold 8 Fold 9 Fold 10 Fold 1 Fold 2 Fold 3 Fold 4 Fold 5 Fold 6 Fold 7 Fold 8 Fold 9 Fold 10 Figura 2.3: Cada recta de regresi´on est´a calculada a partir del subconjunto de entrenamiento de la base de datos que asignamos con la CV . Realizamos la misma operaci´on para los restantes 7 modelos, y colocamos en una matriz la suma de los errores de predicci´on al cuadrado, que es lo que nos interesa: 28 > nf <- 8 > nc <- 1 > resCV <- matrix(nrow=nf, ncol=nc, byrow=TRUE) > rownames(resCV)<-c("mod1","mod2","mod3", "mod4", + "mod5", "mod6", "mod7", "mod8") > colnames(resCV)<-c("ms") > resCV[,1] <- c(attr(RES1, "ms"), + attr(RES2, "ms"), + attr(RES3, "ms"), + attr(RES4, "ms"), + attr(RES5, "ms"), + attr(RES6, "ms"), + attr(RES7, "ms"), + attr(RES8, "ms")) > resCV ms mod1 0.617 mod2 0.615 mod3 0.624 mod4 0.628 mod5 0.620 mod6 0.627 mod7 0.634 mod8 0.632 Por lo que el mejor modelo seg´un la validaci´on cruzada realizada es: modelo 2, price ∼kilometer +age +extra1. 2.4. Criterio de Informaci´on de Akaike El criterio de informaci´on de Akaike, AIC, es un criterio relacionado con el criterio Cpde Mallows, aunque m´as general. La idea principal del AIC es maximizar la log-verosimilitud esperada de un modelo determinado, a trav´es del EMV . El hecho de que se denomine criterio de informaci´on es porque est´a ´ıntimamente relacionado con la llamada informaci´on de Kullback - Leibler. Este criterio no busca encontrar el mejor modelo, sino encontrar el modelo, de entre los que compiten, que mejor se ajuste a los datos con los que trabajamos. Est´a definido por: AIC(M) = −2ln L(ˆ βM,ˆσ2) + 2(|M|+ 1) (2.1) 29 donde ln L(ˆ βM,ˆσ2) es el m´aximo valor del logaritmo de la funci´on verosimilitud evaluado en ˆ βM, donde ˆ βMes el EMV del modelo M, ˆσ2=P(yi−ˆ yi)2 n es el EMV de σ2, y (|M|+ 1) el n´umero total de par´ametros en el modelo incluyendo σ2. El primer t´ermino de la expresi´on 2.1 es una medida de bondad de ajuste, pues disminuye al aumentar ˆ βM, y el segundo t´ermino es una penalizaci´on por el n´umero de par´ametros, exactamente igual que los t´erminos del Cpde Mallows. La generalidad que tiene el criterio AIC, viene del hecho de que podemos calcularlo siempre que tengamos una funci´on de verosimilitud. Entre los modelos que compiten, es mejor el que tenga el menor AIC. Proposici´on 2.4.1. Si se tiene normalidad y σ2es conocida entonces el criterio Cpes equivalente al criterio AIC. Demostraci´on. Comencemos reescribiendo la expresi´on del AIC dada en 2.1: ln L(ˆ yi,ˆ β,ˆσ2) = −n 2ln(2π)−n 2ln(ˆσ2)−n 2P(yi−ˆ yi)2(yi−Xβ)0(yi−Xβ) =−n 2ln(2π)−n 2ln(ˆσ2)−n 2 y excluyendo t´erminos que no dependen del n´umero de par´ametros pllegamos a: −2ln L+ 2(|M|+ 1) = n ln(ˆσ2 p)+2p con p= (|M|+ 1). Ahora, suponiendo conocida σ2, observamos que minimizar el AIC es an´alogo a minimizar: n ln ˆσ2 p σ2+ 2p que puede escribirse, suponiendo normalidad (ˆσ2 p≃σ2) n ln 1 + ˆσ2 p−σ2 σ2+ 2p≃ |{z} ln(1+x)≃xpara x peque˜nos nˆσ2 p σ2−n+ 2p. Luego, sustituyendo σ2por un estimador, qued´andonos con los t´erminos que dependen del n´umero de par´ametros pe ignorando las constantes aditivas en las que no aparece p, pues s´olo dependen de n, llegamos a: AIC ≃Pe2 (p) ˆσ2+ 2p=Cp 30 donde Pe2 (p)es la suma de cuadrados de los residuos del modelo con p par´ametros.  Como consecuencia de la Proposici´on 2.4.1, se obtienen expresiones m´as sencillas del AIC. Corolario 2.4.1. En el modelo lineal con errores gaussianos: −2ln L(ˆ βM,ˆσ2) = nlog(ˆσ2) + 1 ˆσ2(y−XMˆ βM)0(y−XMˆ βM) =nlog(ˆσ2) + nˆσ2 ˆσ2=nlog(ˆσ2) + n donde ˆσ2es el EMV de σ2. Podemos ignorar la constante n, pues no depende de p, y escribir: AIC =nlog(ˆσ2) + 2(|M|+ 1). donde ˆσ2no es el estimador insesgado de σ2. Nota 2.4.1. Las principales ventajas que tiene el criterio AIC, por lo cual es tan pr´actico son: No requiere de ninguna tabla para ver el correspondiente valor. Tiene una f´acil implementaci´on. No necesita un nivel de significaci´on arbitrario para elegir entre dos modelos. Nota 2.4.2. El t´ermino de penalizaci´on no depende del tama˜no de la muestra, es decir, que el n´umero de par´ametros que seleccionamos con este criterio, es el mismo tanto para una muestra peque˜na como para una grande. Lo que hace que AIC no sea consistente, es decir, que no se aproxima al modelo correcto conforme aumenta la muestra, como cabr´ıa esperar. 2.4.1. Aplicaci´on Partimos de los mismos modelos expuestos en el Ejemplo 2.1.1. A diferencia de otros m´etodos, para usar el AIC no necesitamos ninguna librer´ıa extra de R, puesto que ya viene implementado en la funci´on AIC. Lo ´unico que tendremos que hacer es cargar la base de datos y los modelos: 31 > golf <- read.table("golffull.txt", header=TRUE) > attach(golf) > mod1 <- lm(price~kilometer+age, data=golf) > mod2 <- lm(price~kilometer+age+extras1, data=golf) > mod3 <- lm(price~kilometer+age+extras2, data=golf) > mod4 <- lm(price~kilometer+age+TIA, data=golf) > mod5 <- lm(price~kilometer+age+extras1+extras2, data=golf) > mod6 <- lm(price~kilometer+age+extras1+TIA, data=golf) > mod7 <- lm(price~kilometer+age+extras2+TIA, data=golf) > mod8 <- lm(price~kilometer+age+extras1+extras2+TIA, data=golf) Y ahora, basta con calcular los AIC de los modelos y compararlos para ver cual es el mayor. Colocaremos los resultados en una matriz como en ejemplos anteriores para verlo mejor: > nf <- 8 > nc <- 1 > resAIC <- matrix(nrow=nf, ncol=nc, byrow=TRUE) > rownames(resAIC)<-c("mod1","mod2","mod3", "mod4", + "mod5", "mod6", "mod7", "mod8") > colnames(resAIC)<-c("AIC") > resAIC[,1] <- c(AIC(mod1), + AIC(mod2), + AIC(mod3), + AIC(mod4), + AIC(mod5), + AIC(mod6), + AIC(mod7), + AIC(mod8)) > resAIC AIC mod1 406.5904 mod2 405.3859 mod3 408.5433 mod4 408.4380 mod5 407.3697 mod6 406.9542 mod7 410.4026 mod8 408.9490 Por lo que el mejor modelo seg´un este criterio ser´ıa el modelo 2, price ∼ kilometer +age +extra1. 32 2.5. Criterio de Informaci´on Bayesiana Schwarz (1978) ide´o el criterio de informaci´on bayesiana, BIC, a ra´ız de la inconsistencia del estimador AIC, en ´este se considera el tama˜no de la muestra nen el t´ermino de penalizaci´on. Por lo tanto, fue dise˜nado con el objetivo de ser consistente, es decir, que a medida que el tama˜no muestral aumenta, el criterio tiende a seleccionar el verdadero modelo que genera los datos. BIC se basa en estudiar el comportamiento de la probabilidad a posteriori del modelo j-´esimo. Suponiendo ciertas hip´otesis de las distribuciones a priori de los par´ametros, tenemos: ln f(X|Mj) = Lj(ˆ βj|X) + ln P(ˆ βj|Mj) + pj 2ln(2π)−pj 2ln(n) + 1 2ln|Rj| donde f(X|Mj) es la verosimilitud marginal de los datos en el modelo Mj, Lj(ˆ βj|X) es la funci´on de verosimilitud del modelo Mj, que tiene como EMV de βjaˆ βj,P(ˆ βj|Mj) es la probabilidad a priori de los par´ametros, pjes el n´umero de parametros estimados y Rjes igual a nSj, siendo Sjla matriz de covarianzas de ˆ βj. Haciendo tender na infinito, se puede aproximar ln f(X|Mj)≃Lj(ˆ βj|X)−pj 2ln(n) que es equivalente a: BIC(Mj) = −2ln Lj(ˆ βj|X) + ln(n)(pj). Con lo que finalmente llegamos a la expresi´on siguiente: BIC(M) = −2ln L(ˆ β|M|+1) + ln(n)(|M|+ 1).(2.2) Seg´un este criterio, se obtiene el mejor modelo calculando los distintos BIC y qued´andonos con el que tenga el menor. Notemos que si a˜nadimos m´as par´ametros en el modelo, el ajuste se ver´a incrementado, pues el primer t´ermino mide la desviaci´on del modelo estimado con el modelo saturado (con todas las variables), pero este efecto se compensa con el segundo t´ermino pera evitar el sobreajuste. Corolario 2.5.1. Bajo el supuesto de errores gaussianos, la expresi´on 2.2 se reduce a: BIC(M) = nlog(ˆσ2) + log(n)(|M|+ 1). 33 Consideraremos un dato como at´ıpico cuando no se genere por el mismo procedimiento que el resto de la muestra. Por ejemplo, cuando haya un error de medida o si esa observaci´on tiene un valor diferente del resto para una variable explicativa relevante omitida en el modelo. En ese caso, el modelo para esa observaci´on ser´ıa: yi=x0 iβ+ω+ui, donde ωes el error de medida o el efecto de la variable explicativa omitida. Podemos modelizar el dato at´ıpico como un desplazamiento en la media de la distribuci´on. O alternativamente, como un desplazamiento en la varianza, de manera que la observaci´on se genera con nuestro modelo, pero la varianza en ese punto ser´a c2σ2con c >> 1. Ambos modelos son equivalentes, pues con un s´olo dato no es posible saber si la media o la varianza ha cambiado. Un dato at´ıpico puede o no ser influyente, y viceversa. Para buscar los valores at´ıpicos se calculan los residuos estudentizados, ˆ tj, en todos los puntos. Pues el estad´ıstico tasociado a ˆωes precisamente ˆ tj. Definici´on 2.7.2. Se define el residuo estudentizado como: ˆ ti=ei ˆsri√1−hii Para ver si existen valores at´ıpicos se toma el m´aximo de los residuos estudentizados con H0:todos los datos han sido generados por el mismo modelo, ´este seguir´a la distribuci´on del m´aximo de una variable tde Student, que depende de los grados de libertad de ty est´a tabulada. Otro m´etodo para ver ´esto es el m´etodo de Bonferroni y utilizar contrastes m´ultiples. Dicho procedimiento es simple y general, pero no es siempre ´optimo. Se basa en la desigualdad de Bonferroni. Y se utiliza de la forma siguiente: Sea cel n´umero de comparaciones que construimos, sea Aiel suceso: aceptar µi6=µjcuando realmente µi=µj. Supongamos que hacemos las comparaciones de medias con un nivel de significaci´on α: P(Ai) = α. Sea B=A1+A2+... +Ac. Los sucesos Aino son mutuamente excluyentes, por tanto: P(B) = P(A1+A2+... +Ac)≤XP(Ai) = cα. El m´etodo pretende garantizar un error de tipo Itotal para el conjunto de contrastes, αT, por lo que P(B)≤αT.´ Esto se consigue calculando cada contraste individual a un nivel αde manera que: α=αT c. 40 Lo que nos lleva a un procedimiento de aproximaci´on bastante ´util en la pr´actica. Cuando ces grande, se necesitan niveles de significaci´on muy peque˜nos, tanto que no est´an tabulados, por lo que se utiliza la aproximaci´on: tα ν≃1−zα+ 1 4ν−1 donde νson los grados de libertad de tyzαel valor de la distribuci´on normal est´andar (0,1) tal que P(z≥zα) = α. Aunque si tenemos en nuestra muestra un grupo de datos at´ıpicos pueden no detectarse con los procemientos vistos hasta ahora. A pesar de eliminar uno de los puntos, al haber otros parecidos hace que el punto eliminado no parezca influyente. Este fen´omeno se llama enmascaramiento y se resuelve con la estimaci´on robusta y t´ecnicas m´as avanzadas. 2.7.4. Heterocedasticidad Decimos que existe heterocedasticidad en las perturbaciones uicuando no se puede aplicar la hip´otesis: V ar(ui) = σ2, i = 1, .., n con lo que inclumplimos una de las hip´otesis b´asicas donde se asienta la regresi´on lineal. En este caso las observaciones con la varianza baja son importantes, pues son m´as fiables a la hora de estimar la recta de regresi´on que las observaciones con varianza alta (en general, cuanto menor es su varianza, menos se desv´ıan del valor medio que queremos estimar), y deber´ıan tener m´as peso. Pero el m´etodo de m´ınimos cuadrados no tiene en cuenta esto, por lo que los estimadores calculados con este procedimiento dejan de ser eficientes y las f´ormulas deducidas para calcular las varianzas de los estimadores ya no son correctas, por lo tanto, los contrastes basados en ellas dejan de ser v´alidos. La p´erdida de eficiencia de los estimadores depende de la magnitud de heterocedasticidad. Podemos medirla calculando el cociente entre la varianza m´axima y la m´ınima de las observaciones, Bloch y Moses, (1988) recomiendan que cuando el cociente es menor que dos, podemos seguir utiliz´andolos puesto que la p´erdida de eficiencia es peque˜na. Cuando es mayor que dos, la p´erdida de eficiencia es grande. Si adem´as de heterocedasticidad tenemos observaciones con alto efecto palanca, las consecuencias se agravan, pues tambi´en es m´as complicado estimar las perturbaciones del modelo, con lo que es m´as dif´ıcil estimar la varianza 41 de la muestra. Para reconocer la heterocedasticidad basta con analizar los residuos. Mediante el gr´afico de ei=f(ˆyi) se puede detectar, y para identificar si la heterogeneidad en la variabilidad es debida a alguna variable explicativa podemos usar ei=f(xi). Uno de los contrastes para la heterocedasticidad es el de la raz´on de verosimilitudes. Para aplicar este contraste, dividimos los residuos, eien ggrupos, cada uno de un tama˜no niy estimamos la varianza en cada uno de ellos. Sea ˆσ2 ila estimaci´on de la varianza del grupo i, y σ2 iel EMV de la varianza de los residuos. Entonces tenemos el contraste: H0:ei∼N(0, σ) H1:ei∼N(0, σi),con gvalores distintos de σi luego el logaritmo de la raz´on de verosimilitudes de ambas hip´otesis es: log(λ) = − g X i=1 ni 2log(ˆσ2 i)− g X i=1 ni 2−−n 2log(ˆσ2)−n 2 por tanto, 2log(λ) = nlog(ˆσ2)− g X i=1 nilog(ˆσ2 i) cuya distribuci´on asint´otica es χ2 g−1. El contraste anterior no tiene en cuenta la posibilidad de que los residuos sean sesgados por la heterocedasticidad. Para realizar un contraste m´as exacto en muestras peque˜nas, tenemos el siguiente test: H0:yi=x0 iβ+ui, ui∼N(0, σ) H1:yi=x0 iβ+ui, ui∼N(0, σi) que hace que las regresiones sean calculadas por separado en cada grupo al estimar las varianzas, ˆσ2 i. Una vez definido este contraste, se procede de manera an´aloga al anterior. El problema m´as b´asico que produce la heterocedasticidad es la formulaci´on err´onea del modelo. Por ejemplo, si nuestro modelo fuera: y=kxα1 1·... ·xαk k·u donde usigue una distribuci´on log-normal de media 1 y varianza desconocida. Estimamos por un modelo lineal ˆy=ˆ β0+ˆ β1x1+...+ˆ βkxk, los residuos tendr´an falta de normalidad, falta de linealidad y heterocedasticidad, aumentando la varianza de los errores conforme aumentan los valores de las variables explicativas. En este caso, deber´ıamos transformar la variable objetivo, y, con logaritmos. La heterocedasticidad m´as frecuente es que varianza aumente linealmente con el valor de y. Aqu´ı tambi´en se resuelve usando los logaritmos. 42 Si estamos en el caso donde la heterocedasticidad viene por una variable explicativa, xk, y la desviaci´on t´ıpica aumenta linealmente con dicha variable, el procedimiento a seguir es ajustar el siguiente modelo: ˆy xk =ˆ β0 xk +ˆ β1 x1 xk +ˆ βk+u xk donde la perturbaci´on ahora s´ı tiene varianza constante. Otra herramienta ´util para solucionar los problemas de heterocedasticidad es la de m´ınimos cuadrados generalizados. Partamos de un modelo con heterocedasticidad en el que suponemos que: E[UU0] = σ2G donde Ges una matriz sim´etrica y definida positiva. En el caso que nos ocupa, para que las perturbaciones sean heteroced´asticas, se supone que G es una matriz diagonal. Ahora tenemos que diferenciar dos casos, cuando Gsea conocida y cuando no. 1.- Si Ges conocida entonces, Y∼Nn(Xβ, σ2G) y podremos estimar los par´ametros por el m´etodo de m´axima verosimilitud. Que es equivalente a hacer una transformaci´on de las variables con el fin de que cumplan las hip´otesis del modelo de regresi´on y luego aplicar los resultados ya dados. Como Gse supone conocida y definida positiva, podemos obtener una matriz sim´etrica, no singular, Atal que G=AA. Esta Ase denomina matriz ra´ız cuadrada de Gy en nuestro caso, su diagonal son los t´erminos σi/σ. Multiplicando por la inversa de Anuestro sistema, tenemos que: A−1Y=A−1Xβ +A−1U. Esta expresi´on puede reescribirse como: Y∗=X∗β+U∗ con Y∗=A−1Y,X∗=A−1XyU∗=A−1U. Observamos que esta nuevas variables est´an relacionadas entre ellas con el mismo β. Luego la nueva matriz de covarianzas es: E[U∗U∗0] = A−1E[UU0]A−1=σ2I. Con lo que queda arreglado el problema de la heterocedasticidad, pues el modelo Y∗=X∗β+U∗, es homoced´astico, y ya podr´ıamos aplicar el m´etodo 43 de m´ınimos cuadrados (que coincide con el de m´axima verosimilitud) para calcular un estimador de β, que ser´a: ˆ βG= (X∗0X∗)−1X∗0Y∗= (X0G−1X)−1X0G−1Y | {z } G−1=A−1A−1 y se denomina estimador de m´ınimos cuadrados generalizados oMCG y tiene como matriz de covarianzas: V ar(ˆ βG) = σ2(X∗0X∗)−1=σ2(X0G−1X)−1. 2.- Veamos ahora el caso donde Ges desconocida. En general, no es posible resolver el caso en que todos los valores de la matriz Gson desconocidos. Lo habitual es suponer alguna estructura para G, modelizar esta matriz introduciendo par´ametros desconocidos adicionales de forma que el problema planteado sea tratable y utilizar m´etodos iterativos de estimaci´on para el vector de par´ametros βy los nuevos par´ametros utilizados para modelizar la estructura de G. Detalles adicionales pueden verse en Pe˜na, D. [2] (Cap. 9). 2.7.5. Multicolinealidad La multicolinealidad se da cuando las variables explicativas tienen una dependencia entre ellas fuerte, por tanto, es muy dif´ıcil ver el efecto que tiene cada una individualmente en la variable respuesta. Este problema viene del hecho de intentar extraer m´as informaci´on de los datos que lo que contienen, por lo que dicho problema reside en la base de datos y no en el modelo. En los modelos de regresi´on m´ultiple, para estimar el efecto de una variable explicativa debemos fijarnos en la parte de la variable que no est´a relacionada linealmente con las dem´as del modelo. En el caso de que s´ı lo estuviera no ser´ıa posible estimar su efecto, a esto se le llama el problema de la multicolinealidad. Cuando nos disponemos a estimar los par´ametros de los modelos de regresi´on, es necesario invertir la matriz X0X. Si tenemos una variable linealmente dependiente con el resto, la matriz Xtendr´a un rango menor que k+1, que es el n´umero de par´ametros, el determinante de X0Xser´a 0, por lo que la matriz no tendr´a inversa y el sistema de ecuaciones determinado por los par´ametros del modelo tendr´a infinitas soluciones. Puede darse tambi´en que las variables est´en altamente correladas, sin ser exactamente combinaci´on lineal de ninguna, en ese caso habr´a una multicolinealidad alta, por ejemplo, en el caso de que tuvieramos dos variables 44 explicativas en nuestro modelo, x1, x2con medias nulas, tal que: X0X=Px2 1Px1x2 Px1x2Px2 2=ns2 1s12 s12 s2 2 invirtiendo la matriz y utilizando que s12 =rs1s2y|X0X|=s2 1s2 2(1 −r2), tenemos: (X0X)−1=1 n"1 s2 1(1−r2) −r s1s2(1−r2) −r s1s2(1−r2) 1 s2 2(1−r2)#. Luego las varianzas de los estimadores ser´an: V ar(ˆ βi) = σ2 ns2 i(1 −r2), i = 1,2 por tanto, cuando r2∼1 la varianza de los coeficientes estimados ser´a muy alta. Adem´as, las estimaciones tendr´an una gran dependencia entre ellas, pues: Cov(ˆ β1,ˆ β2) = −rσ2 ns1s2(1 −r2). El coeficiente de correlacion entre ˆ β1yˆ β2ser´a igual en valor absoluto, pero de signo contrario, a la correlaci´on entre las variables explicativas, es decir r(ˆ β1,ˆ β2) = Cov(ˆ β1,ˆ β2) qV ar(ˆ β1)qV ar(ˆ β2) =−r. Luego, las estimaciones ser´an tan dependientes entre s´ı, como lo sean las variables entre ellas. En general, la varianza de un coeficiente de regresi´on es V ar(ˆ βi) = σ2/SCR(xi,R) siendo SCR(xi,R) = Pn j=1(xij −ˆxij,R)2la varianza residual de una regresi´on de xisobre el resto. Se tiene tambi´en que SCR(xi,R) = SCT(xi)−SCE(xi,R) = ns2 i(1 −R2 i,R). Llamando Ri,R al coeficiente de correlaci´on m´ultiple en la regresi´on de xien funci´on del resto de variables, tenemos: V ar(ˆ βi) = σ2 ns2 i(1 −R2 i,R) por lo que si el cuadrado del coeficiente de correlaci´on es cercano a 1, la varianza ser´a muy grande. Para averiguar si tenemos o no multicolinealidad debemos examinar: 45 La matriz de correlaci´on entre las variables explicativas, R, y R−1. Los factores de inflaci´on de la varianza. Las ra´ıces y vectores caracter´ısticos de las matrices X0X, o R. Si tenemos una correlaci´on alta entre variables explicativas es una clara se˜nal de multicolinealidad. Puede ser que haya una relaci´on perfecta entre una de las variables explicativas y el resto y, sin embargo, sus coeficientes de correlaci´on sean bajos. Por ejemplo, supongamos las variables explicativas: x1, ..., xkcon media cero, varianza uno y ortogonales. Y definamos una nueva variable que sea la media de las anteriores, xk+1 = (x1, ..., xk)/k. Luego, V ar(xk+1)=1/k yCov(xi, xk+1)=1/k, con lo que su correlaci´on es 1/√k. Si kes grande, la correlaci´on ser´a peque˜na, pero un modelo que incluya las k+ 1 variables (x1, ..., xk, xk+1) tendr´a una multicolinealidad exacta. Sea Rla matriz de correlaci´on de las variables explicativas, la cual es cuadrada, sim´etrica de orden ky cuyo t´ermino (ij) es el coeficiente de correlaci´on lineal simple entre xiyxj. Esta matriz, para dos variables ser´ıa: R=1r r1 luego R−1es: R−1=1 1−r2 −r 1−r2 −r 1−r2 1 1−r2 podemos ver que los elementos de la diagonal, 1/(1 −r2), contienen al coeficiente de correlaci´on. Para kvariables, los elementos de la diagonal ser´ıan 1/(1 −R2 i,R), siendo R2 i,R el coeficiente de correlaci´on m´ultiple de la variable explicativa xicon el resto de variables explicativas. Por tanto, si tenemos elementos de la diagonal de R−1grandes, nos indicar´a que hay alta multicolinealidad. En este caso no tenemos el problema que ten´ıamos con los elementos de R, donde pod´ıa no detectarse a simple vista la multicolinealidad, pues en los elementos de la diagonal de la matriz inversa se tienen en cuenta todas las variables explicativas, y entonces se detectar´a la multicolinealidad cuando una de las variables sea casi combinaci´on lineal del resto. Aunque R−1tambi´en tiene inconvenientes, cuando la matriz Rsea casi singular, no podremos calcular su inversa con precisi´on. Los t´erminos de la diagonal de R−1se interpretan como el aumento de la variabilidad en la estimaci´on de los efectos de cada variable explicativa en la regresi´on m´ultiple, como consecuencia de la dependencia entre las variables, respecto a la regresi´on simple. Ve´amoslo para dos variables explicativas de media cero: 46 La varianza de las estimaciones de los efectos de las variables mediante regresiones simples ser´ıa ˆs2 r(i)/s2 in, con s2 r(i) la varianza residual de la regresi´on simple que tiene por regresor la variable xi. Si estimamos los efectos mediante regresi´on m´ultiple, la varianza ser´ıa ˆs2 r/s2 i(1 −r2)n. Por tanto: V ar(efecto xi|R. m´ultiple) V ar(efecto xi|R. simple) =ˆs2 r(i) ˆs2 r 1 1−r2 dicha expresi´on nos indica que el cambio de la varianza de un coeficiente al pasar de la regresi´on simple a la regresi´on m´ultiple depende de dos factores. Uno, el cambio de la varianza residual de la regresi´on, que ser´a mayor en la simple que en la m´ultiple, normalmente. Y dos, el 1/(1 −r2), denominado factor de inflaci´on de la varianza, el cual mide el aumento de la varianza debido a la dependencia entre las variables. Por ejemplo, al introducir una nueva variable explicativa en la regresi´on simple la cual est´e muy correlada con la que ya hab´ıa y no ayuda a explicar la variable objetivo, hace que el primer t´ermino, ˆs2 r(i)/ˆs2 r, est´e cercano a 1 y la varianza del coeficiente de la primera variable estar´a multiplicada por el mencionado factor de inflaci´on. Se puede probar que el resultado anterior se puede generalizar como sigue: V ar(efecto xi|R. m´ultiple) V ar(efecto xi|R. simple) =ˆs2 r(i) ˆs2 r FIV (i) donde FIV (i) = 1/(1 −R2 i,R) es el factor de inflaci´on de la varianza. Cuando X0XoRson singulares debemos recurrir a otras t´ecnicas para tratar la multicolinealidad. Por ejemplo: el ´ındice de condicionamiento, denotado por IC, el cual nos sirve para estos casos y est´a definido por: IC =rm´aximo autovalor de la matriz m´ınimo autovalor de la matriz ≥1. Normalmente se calcula este ´ındice para Ren vez de para X0X, puesto que ´esta no est´a afectada por las escalas de los regresores, pues en el caso de que un regresor tuviera una varianza grande y otro muy peque˜na, la matriz X0X estar´ıa mal condicionada, y por tanto, fuera de la diagonal tendr´ıa t´erminos nulos. Por convenio se admite que existe alta multicolinealidad cuando IC > 30. Cuando 10 < IC < 30 tendremos una multicolinealidad moderada. Y en caso contrario tendremos bien definida la matriz y la multicolinealidad ser´a lo suficientemente baja para no alterar la estimaci´on por el m´etodo de m´ınimos cuadrados del modelo. 47 Antes de pasar a ver como solucionar la multicolinealidad, veamos el efecto de ´esta en el error cu´adratico medio, lo que nos ser´a ´util para ver una de sus soluciones. Tenemos que el error cu´adratico medio de ˆ βest´a definido por: ECM(ˆ β) = E[(ˆ β−β)0(ˆ β−β)] = k X i=0 (ˆ βi−βi)2=σ2tr(X0X)−1=σ2 k X i=0 1 λi donde los λison los valores propios de la matriz X0X. Si esta matriz es casi singular, λi≃0 para alg´un i, lo que lleva a tener un error cuadr´atico medio muy grande. Una vez visto esto, pasemos a mostrar como solucionar la multicolinealidad. ´ Esta no tiene soluci´on sencilla pues, como mencionamos al principio de la secci´on, el problema reside en la muestra. Una de la alternativas es tomar las observaciones de manera que la matriz X0Xsea diagonal, lo que reduce la varianza de los estimadores. En caso de no poder dise˜nar la manera de recabar los datos, podemos eliminar regresores altamente correlados con otros, haciendo menor el n´umero de par´ametros a estimar, aunque dichos estimadores ser´an sesgados. Esta es una de las soluciones m´as simples. Veamos la manera de proceder. Sea y=β1x1+β2x2+u nuestro modelo, donde vamos a suponer que las variables explicativas tienen media cero, por simplificar. Seg´un lo visto anteriormente, V ar(ˆ β1) = σ2/ns2 1(1 −r2 12). Por lo que su error cuadr´atico medio es: ECM(ˆ β1) = V ar(ˆ β1). Si eliminamos la variable explicativa x2, nos queda el modelo: y=b1x1+ε por lo que la estimaci´on de b1ser´a: ˆ b1=Pyx1 Px2 1 veamos que en efecto, es sesgada: E[ˆ b1] = 1 Px2 1 Ehβ1Xx2 1+β2Xx2x1+uXx1i 48 =β1+β2Px2x1 Px2 1 =β1+β2r12 s2 s1 por lo que, s´ı, es sesgado. Calculando su varianza, obtenemos: V ar(ˆ b1) = σ2 Px2 1 =σ2 ns2 1 . Y por tanto, su error cuadr´atico medio ser´a: ECM(ˆ b1) = β2r12 s2 s12 +σ2 ns2 1 por lo que deber´ıa verificarse que ECM(ˆ b1)< ECM(ˆ β1), es decir: β2 2r2 12ns2 2+σ2 ns2 1 <σ2 ns2 1(1 −r2 12) de lo que se deduce que, 1 1−r2 12 >β2 σ2 ns2 2 luego cuando r12 ≃1, el ECM(ˆ b1) ser´a menor que ECM(ˆ β1) y obtendremos una estimaci´on mejor (aunque sesgada) del efecto de la variable explicativa x1eliminando de nuestro modelo la variable explicativa x2. Nota 2.7.1. Reordenando la ´ultima expresi´on obtenemos un resultado bastante interesante: 1>β2 2ns2 2(1 −r2 12) σ2=β2 2 V ar(ˆ β2)= β2 DT(ˆ β2)!2 donde DT(ˆ β2)es la desviaci´on t´ıpica de ˆ β2. Sustituyendo los par´ametros β2 yσ2por sus estimaciones, obtenemos el estad´ıstico t, al cuadrado, que se utiliza para contrastar si el par´ametro es cero. Teniendo en cuenta la Nota 2.7.1, eliminaremos de nuestro modelo las variables cuyo estad´ıstico tsea menor que 1, para as´ı tratar de mejorar el error cuadr´atico medio de estimaci´on de los par´ametros restantes y eliminar la multicolinealidad. En vez de eliminar directamente las variables de nuestro modelo, podemos crear una nueva variable que agrupe las que est´an muy correladas entre s´ı. 49 Step: AIC=-84.73 price ~ kilometer + age + extras1 Df Sum of Sq RSS AIC <none> 100.32 -84.729 - extras1 1 1.887 102.21 -83.524 - kilometer 1 27.768 128.09 -44.700 - age 1 76.600 176.92 10.852 Por tanto, seg´un este m´etodo obtenemos el siguiente modelo: price ∼kilometer+ age +extras1. 3.3. Selecci´on paso a paso El m´etodo de selecci´on paso a paso (o stepwise selection), es una combinaci´on de los dos anteriores, as´ı, evita los inconvenientes de la selecci´on hacia adelante y no requiere de una capacidad de c´alculo tan grande como la de la selecci´on hacia atr´as. En cada paso se contrasta si entra una nueva variable explicativa o sale una que ya est´e en el modelo. El algoritmo requiere fijar dos reglas, una para las variables de entrada y otra para las variables de salida. El proceso termina cuando no haya mejoras significativas a la hora de a˜nadir o eliminar alguna variable. Este m´etodo es el m´as utilizado de los 3. 3.3.1. Aplicaci´on Por ´ultimo veremos un ejemplo de la selecci´on paso a paso, en el cual utilizaremos la misma librer´ıa MASS con su funci´on stepAIC, que recibir´a de argumentos: el modelo con el que comienza el algoritmo, el modelo con el m´aximo n´umero de variables, y la direcci´on en la que avanza, en este caso, hacia ambos lados. Como en el ejemplo anterior, ejecutamos solamente la funci´on: > mod.step <- stepAIC(mod0, scope = list(upper = mod8), > direction = "both") Start: AIC=76.69 price ~ 1 Df Sum of Sq RSS AIC + age 1 135.435 130.09 -44.030 + kilometer 1 88.086 177.44 9.357 56 + extras2 1 3.663 261.86 76.297 <none> 265.53 76.686 + TIA 1 0.286 265.24 78.501 + extras1 1 0.257 265.27 78.520 Step: AIC=-44.03 price ~ age Df Sum of Sq RSS AIC + kilometer 1 27.886 102.21 -83.524 + extras1 1 2.004 128.09 -44.700 <none> 130.09 -44.030 + extras2 1 0.543 129.55 -42.750 + TIA 1 0.200 129.89 -42.294 - age 1 135.435 265.53 76.686 Step: AIC=-83.52 price ~ age + kilometer Df Sum of Sq RSS AIC + extras1 1 1.887 100.32 -84.729 <none> 102.21 -83.524 + TIA 1 0.091 102.12 -81.677 + extras2 1 0.028 102.18 -81.572 - kilometer 1 27.886 130.09 -44.030 - age 1 75.234 177.44 9.357 Step: AIC=-84.73 price ~ age + kilometer + extras1 Df Sum of Sq RSS AIC <none> 100.32 -84.729 - extras1 1 1.887 102.21 -83.524 + TIA 1 0.251 100.07 -83.161 + extras2 1 0.009 100.31 -82.745 - kilometer 1 27.768 128.09 -44.700 - age 1 76.600 176.92 10.852 Con este m´etodo obtenemos el siguiente modelo: price ∼age +kilometer + extras1. 57 58 Cap´ıtulo 4 T´ecnicas de regularizaci´on Para hallar los estimadores por el m´etodo de m´ınimos cuadrados ordinarios de los par´ametros en el modelo lineal cl´asico, hay que resolver el sistema de ecuaciones: X0Xβ =X0y. Para que este sistema tenga soluci´on ´unica, la matriz Xdebe tener rango m´aximo, rango(X) = p. Sin embargo, puede haber problemas: Cuando haya columnas en la matriz Xque sean casi combinaci´on lineal de otras, es decir, cuando se presenta el problema de la colinealidad. Cuando el n´umero de regresores, p, es grande. En este caso, la soluci´on es num´ericamente inestable, aunque los coeficientes sigan siendo identificables en teor´ıa. Adem´as, en muchas de las aplicaciones que est´an surgiendo en nuestros d´ıas, por ejemplo en gen´etica, ocurre que el n´umero de covariables es mucho mayor que el n´umero de observaciones de las que disponemos. Esto se conoce como problemas con “n peque˜no, y p grande”. En todas estas situaciones son ´utiles las t´ecnicas de regularizaci´on. Puede decirse que las t´ecnicas de regularizaci´on se aplican para obtener estimaciones de los coeficientes de regresi´on cuando la matriz X0Xes singular o est´a muy pr´oxima a serlo. En este contexto, regularizar significa, hacer el problema tratable, imponiendo una serie de restricciones al conjunto de soluciones admisibles. En las t´ecnicas de regularizaci´on se plantea un problema de optimizaci´on que considera como funci´on objetivo una obtenida por m´ınimos cuadrados penalizados (Penalized Least Squares, PLS) P LS(β) = (y−Xβ)0(y−Xβ)+λpen(β) 59 donde λ≥0 es un par´ametro de penalizaci´on que controla el efecto de la penalizaci´on y pen(β) es el t´ermino de penalizaci´on. Si λ≃0 entonces ˆ βP LS est´a pr´oximo al MCO( ˆ βLS). En cambio, si λes grande, se le da mucha importancia a la penalizaci´on. Luego, se plantea ahora el problema de hallar ˆ βP LS: ˆ βP LS = arg min β [(y−Xβ)0(y−Xβ)+λpen(β)]. Que es equivalente a resolver el problema de optimizaci´on: ˆ βP LS = arg min β [(y−Xβ)0(y−Xβ)] s.a. pen(β)≤t donde t es una constante relacionada con el par´ametro de penalizaci´on (smoothing) λen una relaci´on uno a uno. 4.1. Regresi´on contra´ıda La regresi´on contra´ıda (o ridge regression) fue introducida en 1970 por Hoerl y Kennard. Recordemos que β= (β1, ..., βp)0. En el caso de la regresi´on ridge se considera la siguiente penalizaci´on: pen(β) = kβk2= k X j=0 β2 j=β0β por tanto nos queda: PLS(β) = (y−Xβ)0(y−Xβ)+λβ0β. De donde puede comprobarse que: ˆ βP LS = (X0X+λIp)−1X0y. Observaci´on: Recordemos que ˆ βLS = (X0X)−1X0y. 60 Para valores de λcercanos a cero, el impacto de pen(β) es pr´acticamente nulo, y ˆ βP LS ≃ˆ βLS. Sin embargo, si λes grande, esta t´ecnica permite resolver el problema de la multicolinealidad, porque hace que la matriz (X0X+λIp)−1 sea invertible en el caso de que X0Xno lo fuera. Adem´as, ˆ βP LS es una contracci´on de ˆ βLS hacia cero. Esto puede verse observando la funci´on objetivo a minimizar, y el papel que en ella desempe˜na λβ0β. Si λes grande, el min β{PLS(β)}, estar´a determinado por el t´ermino que minimice λpen(β) = λβ0β, que claramente se minimiza cuando β= 0. En la pr´actica no interesa penalizar la ordenada en el origen (intercept) del modelo de regresi´on, β0. Para ello existen dos alternativas: Centrar todas las covariables y la variable respuesta, Y, para que los nuevos valores de ´estas tengan media cero, y= 0, x= 0, lo que autom´aticamente produce que ˆ β0= 0. Esto implica que el intercept (o t´ermino constante) se elimina del modelo, y por lo tanto no se penaliza. Modificar la penalizaci´on a: pen(β) = k X j=1 β2 j=β0Kβ donde K=diag(0,1, ..., 1), es decir, se introduce una matriz de penalizaci´on que excluye al coeficiente β0, y sigue siendo la identidad para el resto de los coeficientes. Esta segunda opci´on es la que adoptaremos, y nos conduce al estimador de regresi´on contra´ıda (ridge estimate) dado por: ˆ βP LS = (X0X+λK)−1X0y. Proposici´on 4.1.1. Propiedades de ˆ βP LS: E(ˆ βP LS)=(X0X+λK)−1X0Xβ. Por tanto, no es insesgado salvo que λ= 0. Normalmente ocurrir´a |ˆ βj,P LS| ≤ |ˆ βj,LS|,j= 1, ..., k. Aunque no siempre se tiene. Cov(ˆ βP LS) = σ2(X0X+λK)−1X0X(X0X+λK)−1. Proposici´on 4.1.2. Comparaci´on ˆ βP LS,ˆ βLS: E(ˆ βLS) = E((X0X)−1X0y) = (X0X)−1X0Xβ =β, luego es insesgado. Cov(ˆ βLS) = σ2(X0X)−1. 61 Luego: ˆ βP LS = (X0X+λK)−1X0y= (X0X+λK)−1(X0X) (X0X)−1X0y | {z } ˆ βLS = (X0X+λK)−1(X0X)ˆ βLS. En el caso de matrices ortogonales se puede obtener una relaci´on entre los coeficientes de ambos estimadores que ilustra porqu´e se llama regresi´on contra´ıda: Cov(ˆ βP LS) = σ2(X0X+λK)−1(X0X)(X0X+λK)−1. Puede probrarse que la matriz Cov(ˆ βLS)−Cov(ˆ βP LS) es definida positiva para λ > 0, lo que implica que: V ar(ˆ βj,P LS)< V ar(ˆ βj,LS), j = 1, ..., k. En resumen, con la regresi´on ridge el estimador que se obtiene es sesgado, pero tiene menor ECM. Lo ´unico que quedar´ıa es calcular el par´ametro λadecuado, lo que se hace normalmente por el m´etodo de la validaci´on cruzada comentado en la Secci´on 2.3. N´otese por ´ultimo que la escala de las covariables es importante cuando estamos regularizando. La penalizaci´on formada por el cuadrado de los coeficientes de regresi´on asume que todos los coeficientes pueden ser comparados en valor absoluto. Sin embargo, la escala tiene un impacto directo en la interpretacion de esos valores absolutos. Por ejemplo, el coeficiente asociado a una covariable que mide la distancia estar´a escalada por un factor de 1.000 cuando la variable se mida en metros en vez de kil´ometros. Por lo tanto, es importante hacer que todas las variables sean comparables en su escala antes de aplicar la aproximaci´on por m´ınimos cuadrados penalizados. La soluci´on m´as com´un es la de normalizar todas las variables. 4.1.1. Aplicaci´on Para ilustrar la regresi´on ridge usaremos los paquetes de R:car yMASS. Y un fichero de datos utilizado en Tibshirani, R. et al. [12]. El objetivo de este estudio es determinar qu´e variables influyen en la presencia de un ant´ıgeno prost´atico espec´ıfico, el cual se utiliza para detectar el c´ancer de pr´ostata. > url <- "http://www-stat.stanford.edu/~tibs/ElemStatLearn + /datasets/prostate.data" > cancer <- read.table(url, header=TRUE) > library(car) > library(MASS) 62 Nuestro fichero de datos cuenta con 97 observaciones y 10 variables, las cuales son: lcavol: log-vol´umen del c´ancer. lweight: log-tama˜no de la pr´ostata. age: edad del paciente. lbhp: log-cantidad de hiperplasia benigna. svi: toma el valor 1 si est´a invadida la ves´ıcula seminal y 0 si no. lcp: log-penetraci´on capsular. gleason: puntuaci´on Gleason. pgg45: porcentaje de la puntuaci´on Gleason 4 ´o 5. lpsa: log-an´alisis del ant´ıgeno prost´atico espec´ıfico. train: variable para distinguir el conjunto de entrenamiento y el de test. Seleccionemos ahora el conjunto test y el conjunto train vali´endonos de la variable ”train” antes mencionada. Tal como est´a conformado el fichero de datos, el 70 % del conjunto est´a destinado al entrenamiento del modelo. > train = subset(cancer,train=="TRUE") > test = subset(cancer,train=="FALSE") Calculemos ahora el modelo de regresi´on ridge con la funci´on lm.ridge la cual tiene implementada una b´usqueda del λ´optimo a trav´es de la validaci´on cruzada generalizada, es importante remarcar que este t´ermino puede inducir a error, pues no es una generalizaci´on de la validaci´on cruzada mencionada en la secci´on 2.3, aunque se utiliza por convenio dicho nombre, se podr´ıa hablar de “aproximaci´on”. > modelo_ridge <- lm.ridge(lpsa ~ ., data=train[,-10], > lambda = seq(0,10,0.1)) > plot(seq(0,10,0.1), modelo_contraida$GCV, > main="B´usqueda lambda por GCV", + type="l", xlab=expression(lambda), ylab="GCV") 63 0246810 0.00832 0.00834 0.00836 0.00838 0.00840 0.00842 0.00844 Búsqueda lambda por GCV λ GCV Vemos que el λdebe estar pr´oximo a 5, para averiguar el valor ´optimo podemos emplear la funci´on select: > select(lm.ridge(lpsa ~ ., data=train[,-10], lambda = seq(0,10,0.1))) modified HKB estimator is 3.355691 modified L-W estimator is 3.050708 smallest value of GCV at 4.9 luego el valor ´optimo es λ= 4,9. Podemos ver tambi´en c´omo var´ıan los coeficientes al modificar el λ > matplot(seq(0,10,0.1), coef(modelo_ridge)[,-1], xlim=c(0,11), type="l", + xlab=expression(lambda), ylab=expression(hat(beta)), lty=1, lwd=2, + main="Coeficientes en funci´on del lambda") > text(rep(10, 9), coef(modelo_ridge)[length(seq(0,10,0.1)),-1], + colnames(train)[-9], pos=4) 64 0246810 −0.2 0.0 0.2 0.4 0.6 Coeficientes en función del lambda λ β ^ lcavol lweight age lbph svi lcp gleason pgg45 train Se aprecia que al aumentar el λlos coeficientes tienden a 0, pero debemos tener en cuenta que a mayor λ, mayor es el sesgo de nuestro modelo. Tenemos ya definido nuestro modelo: > modelo_ridge <- lm.ridge(lpsa ~ ., data=train[,-10], lambda = 4.9) > coefficients(modelo_ridge) lcavol lweight age lbph 0.096814771 0.492787412 0.601103227 -0.014821787 0.138019854 svi lcp gleason pgg45 0.679632580 -0.116790333 0.017113954 0.007081258 Para finalizar calculemos el error cuadr´atico medio del modelo obtenido por el m´etodo de m´ınimos cuadrados ordinarios y el error del modelo penalizado. > modelo_mco <- lm(lpsa~ . , data=train[,-10]) > ajuste_mco <- predict(modelo_mco,test) > sum((test$lpsa-ajuste_mco)^2) [1] 15.63822 Este modelo tiene una suma de errores cuadr´aticos medios de 15.63822. Veamos ahora el modelo obtenido por la regresi´on ridge: > coeficientes <- as.vector(coef(modelo_ridge)) > matriz <- as.matrix(test[,-9:-10]) > matriz <- cbind(rep(1,length=nrow(test)),matriz) 65 penalizados. Para ilustrar este hecho, vamos a considerar un vector de coeficientes β= (β1, β2)0; pero todos los resultados se generalizan f´acilmente. Notemos que no hemos incluido al intercept lo que supone que consideremos las covariables estandarizadas y una variable respuesta centrada. Proposici´on 4.3.1. El criterio de m´ınimos cuadrados, LS(β), puede reescribirse como: LS(β)=(β−ˆ β)0X0X(β−ˆ β) + y0(In−X(X0X)−1X0)y =(β−ˆ β)0X0X(β−ˆ β) + ˆε0ˆε.(4.1) Demostraci´on. Partamos del producto del criterio de m´ınimos cuadrados (y−Xβ)0(y−Xβ) = y0y−2β0X0y+β0X0Xβ. Veamos ahora la expansi´on de la forma cuadr´atica en β: (β−ˆ β)0X0X(β−ˆ β) = β0X0Xβ −2β0X0Xˆ β+ˆ β0X0Xˆ β. Teniendo en cuenta que ˆ β= (X0X)−1X0y, el segundo sumando es 2β0X0Xˆ β= 2β0X0X(X0X)−1X0y= 2β0X0y y el tercer sumando es ˆ β0X0Xˆ β=y0X(X0X)−1X0X(X0X)−1X0y=y0X(X0X)−1X0y. Con lo que llegamos a (β−ˆ β)0X0X(β−ˆ β) = β0X0Xβ −2βX0y+y0X(X0X)−1X0y con esto tendr´ıamos el primer sumando, el segundo la obtenemos viendo en la Expresi´on 4.1 que ˆε0ˆε=(y−Xˆ β)0(y−Xˆ β) =(y−X(X0X)−1X0y)0(y−X(X0X)−1X0y) =y0(In−X(X0X)−1X0)y.  Como hemos visto que LS(β) es equivalente a una forma cuadr´atica, los valores de βque resuelven LS(β) = c, para una constante c, es decir, sus curvas de nivel, son elipses con una forma determinada por la matriz X0X. Por otro lado, en dos dimensiones, la restricci´on: |β1|+|β2|=t 72 define curvas de nivel con forma de diamente de lado √2t. Por lo tanto, el estimador LASSO regularizado, dada una t, es el punto de corte de las dos regiones geom´etricas definidas por la restricci´on y por el criterio de los m´ınimos cuadrados. Si el punto de corte est´a en uno de los v´ertices del diamante, algunos coeficientes se estimar´an como cero. Las curvas de nivel que se definen en la regresi´on ridge son c´ırculos de la forma: β2 1+β2 2=t con lo que no se puede dar el caso que el estimador corte con la regi´on en un v´ertice, pues es un c´ırculo, con lo que no conseguiremos que ning´un coeficiente se estime como cero. Lo vemos mejor en los gr´aficos recogidos en la Figura 4.3: 73 -4-2 0 2 4 6 -4 -2 0 2 4 6 beta1 beta2 (a) -4-2 0 2 4 6 -4 -2 0 2 4 6 beta1 beta2 (b) -4-2 0 2 4 6 8 -4 -2 0 2 4 6 8 beta1 beta2 (c) -4-202468 -4 -2 0 2 4 6 8 beta1 beta2 (d) Figura 4.3: Interpretaci´on geom´etrica del criterio de m´ınimos cuadrados penalizados para la regresi´on ridge (figuras (a) y (c)) y para la regresi´on LASSO (figuras (b) y (d)). En las figuras superiores, se considera una matriz X0X no diagonal, mientras que las figuras inferiores corresponden a una matriz de dise˜no X0X=I2. Resumen M´ınimos cuadrados penalizados La estimaci´on regularizada en el modelo lineal permite penalizar el criterio de los m´ınimos cuadrados: PLS(β) = (y−Xβ)0(y−Xβ)+λpen(β) con el par´ametro de penalizaci´on, λ≥0. 74 Regresi´on ridge Para la regresion ridge, la penalizaci´on viene dada por la suma de los coeficientes al cuadrado: pen(β)= k X j=1 β2 j=β0Kβ con la matriz de penalizaci´on K=diag(0,1, ..., 1). El resultado de la estimaci´on de los m´ınimos cuadrados penalizados es: ˆ βP LS = (X0X+λK)−1X0y. Regresi´on LASSO Para la regresi´on LASSO, la penalizaci´on viene dada por la suma de los valores absolutos de los coeficientes: pen(β) = k X j=1 |βj|. El resultado de la estimaci´on no tiene una f´ormula anal´ıtica y debe determinarse num´ericamente, por ejemplo utilizando t´ecnicas de programaci´on cuadr´atica. Elecci´on del par´ametro de penalizaci´on El par´ametro de penalizaci´on, λ, puede determinarse por los m´etodos r-fold cross validation o con la generalized cross validation. 75 Ap´endice A Anexo A.1. Comandos en Rde las gr´aficas A.1.1. Figura 2.1 > par(mfrow=c(2,2), las=1) > plot(age,price,ylab="sales price in 1000$",xlab="age in months", + main="Sales price vs age") > plot(kilometer,price,ylab="sales price in 1000$", + xlab="kilometer reading in 1000 km", main="Sales price vs kilometer") > plot(TIA,price,ylab="sales price in 1000$", + xlab="months until next TIA appointment", main="Sales price vs TIA") A.1.2. Figura 2.2 > par(mfrow=c(1,2), las=2) > boxplot(price ~ extras1, main="Sales price vs no ABS/ABS", + ylab="sales in 1000$", col="gold", names=c("no ABS", "ABS")) > boxplot(price ~ extras2, main="Sales price vs no sunroof/sunroof", + ylab="sales in 1000$", col="gold", names=c("no sunroof", "sunroof")) A.1.3. Figura 4.1 > library(gplots) > curve(x^2, from=-2, to=2, xlab=expression(beta), + ylab=expression(pen(beta)), col="red", ylim=c(0, 4)) > curve(abs(x), from=-2, to=2, col="blue", add=T) > legend("topright", c("Ridge", "Lasso"), + lwd=2, col=c("red", "blue")) 76 A.2. Comandos en Mathematica de las gr´aficas A.2.1. Figura 4.3 Penalizaci´on Ridge y matriz no -diagonal circ1 := ParametricPlot[{r*Cos[t], r*Sin[t]}, {t, 0, 2*Pi}, {r, 0, 2.16}, PlotStyle -> {GrayLevel[0.25]}, PlotRange -> {{-5, 7}, {-5, 7}}, AxesLabel -> {"beta1", "beta2"}] circ2 := ParametricPlot[{r*Cos[t], r*Sin[t]}, {t, 0, 2*Pi}, {r, 2.16, 3.32}, PlotStyle -> {GrayLevel[0.55]}, PlotRange -> {{-5, 7}, {-5, 7}}, AxesLabel -> {"beta1", "beta2"}] circ3 := ParametricPlot[{r*Cos[t], r*Sin[t]}, {t, 0, 2*Pi}, {r, 3.32, 4.84}, PlotStyle -> {GrayLevel[0.8]}, PlotRange -> {{-5, 7}, {-5, 7}}, AxesLabel -> {"beta1", "beta2"}] X = {{0.5, -1}, {0.45, 0.5}}; beta = {{beta1}, {beta2}}; LS = (beta - {{8}, {3}})\[Transpose].X\[Transpose].X.(beta - {{8}, \ {3}}); elipse := ContourPlot[(beta - {{8}, {3}})\[Transpose].X\[Transpose].X.(beta - \ {{8}, {3}}), {beta1, 0, 7}, {beta2, -2, 5}, ContourShading -> False, PlotRange -> {{-5, 7}, {-5, 7}}, AxesLabel -> {"beta1", "beta2"}] cont1 := Show[elipse, circ1, circ2, circ3, Axes -> True] cont1 Penalizaci´on LASSO y matriz no - diagonal rombo1 := ContourPlot[{Abs[beta1] + Abs[beta2] == 1.4}, {beta1, -4, 4}, {beta2, -4, 4}, PlotRange -> {{-5, 7}, {-5, 7}}, AxesLabel -> {"beta1", "beta2"}] rombo2 := ContourPlot[{Abs[beta1] + Abs[beta2] == 2.6}, {beta1, -4, 4}, {beta2, -4, 4}, PlotRange -> {{-5, 7}, {-5, 7}}, AxesLabel -> {"beta1", "beta2"}] rombo3 := ContourPlot[{Abs[beta1] + Abs[beta2] == 4.2}, {beta1, -4.5, 77 4.5}, {beta2, -5, 5}, PlotRange -> {{-5, 7}, {-5, 7}}, AxesLabel -> {"beta1", "beta2"}] cont2 := Show[elipse, rombo1, rombo2, rombo3, Axes -> True] cont2 Penalizaci´on Ridge y matriz ortonormal X = {{1, 0}, {0, 1}}; circulos := ContourPlot[(beta - {{6}, {3}})\[Transpose].X\[Transpose].X.(beta - \ {{6}, {3}}), {beta1, 0, 8}, {beta2, -5, 10}, ContourShading -> False, PlotRange -> {{-5, 9}, {-5, 9}}, AxesLabel -> {"beta1", "beta2"}] cont3 := Show[circulos, circ1, circ2, circ3, Axes -> True] cont3 Penalizaci´on LASSO y matriz ortonormal cont4 := Show[circulos, rombo1, rombo2, rombo3, Axes -> True] cont4 A.3. Paquetes de R C. Agostinelli and U. Lund (2013). R package ’circular’: Circular Statistics (version 0.4-7). URL https://r-forge.r-project.org/projects/ circular/. Claudio Agostinelli and SLATEC Common Mathematical Library (2015). wle: Weighted Likelihood Estimation. R package version 0.9-91. https: //CRAN.R-project.org/package=wle. John H. Maindonald and W. John Braun (2015). DAAG: Data Analysis and Graphics Data and Functions. R package version 1.22. https: //CRAN.R-project.org/package=DAAG. Sarkar, Deepayan (2008) Lattice: Multivariate Data Visualization with R. Springer, New York. ISBN 978-0-387-75968-5. Venables, W. N. & Ripley, B. D. (2002) Modern Applied Statistics with S. Fourth Edition. Springer, New York. ISBN 0-387-95457-0. 78 John Fox and Sanford Weisberg (2011). An R Companion to Applied Regression, Second Edition. Thousand Oaks CA: Sage. URL: http: //socserv.socsci.mcmaster.ca/jfox/Books/Companion. Jelle Goeman, Rosa Meijer and Nimisha Chaturvedi (2014). penalized: L1 (lasso and fused lasso) and L2 (ridge) penalized estimation in GLMs and in the Cox model. R package version 0.9-45. https: //CRAN.R-project.org/package=penalized. 79 80 Bibliograf´ıa [1] Fahrmeir, L.; Kneib, Th.; Lang, S.; Marx, B. Regression: Models, Methods and Applications. New York: Springer. 2013. [2] Pe˜na, D. Regresi´on y dise˜no de experimentos. Madrid: Alianza Editorial, S.A. 2002. [3] Harrell, F. E. Jr. Regression Modeling Strategies. New York: Springer. 2001. [4] Apuntes de Modelos lineales y dise˜no de experimentos. Tercer curso de Grado en Matem´aticas, 2013-14. Universidad de Sevilla. (Profesores D. Juan M. Mu˜noz Pichardo y D. Joaqu´ın Antonio Garc´ıa de las Heras). [5] Apuntes de Inferencia estad´ıstica. Tercer curso de Grado en Matem´aticas, 2014-15. Universidad de Sevilla. (Profesores D. Emilio Carrizosa Priego y D. Joaqu´ın Antonio Garc´ıa de las Heras). [6] Carmen Garc´ıa Olaverri. (1996). Estabilidad de algunos criterios de selecci´on de modelos. Q¨uestii´o, Vol 20, 2 pp. 147-166. [7] Andrew W. Moore. Cross-validation for detecting and preventing overfitting. Apuntes. Carnegie Mellon University. [8] Pilar Cacheiro Mart´ınez. (2011). M´etodos de selecci´on de variables en estudios de asociaci´on gen´etica. Aplicaci´on a un estudio de genes candidatos en Enfermedad de Parkinson. Fin de m´aster. A Coru˜na: Universidad de Santiago de Compostela. [9] D˜na. Mar´ıa Jes´us B´arcena Ru´ız. Universidad del Pa´ıs Vasco. Econom´ıa Aplicada III (Estad´ıstica y Econometr´ıa). http://campusvirtual. ehu.es/open_course_ware/castellano/experimentales/ estadistica/materiales-de-estudio/index.html. [10] D. Juan M. Vilar Fern´andez. Universidad de Santiago de Compostela. http://dm.udc.es/asignaturas/estadistica2/indice_res.html. 81