Estadística y Machine Learning con R
Francisco Parra
15 de Junio de 2017
Extensiones al modelo de regresión lineal
Introducción
Como veíamos en el capitulo anterior, el modelo de regresión lineal requiere que se cumplan las siguientes hipótesis sobre los términos de error:
Media cero: \(E(e_i) = 0\) \(i=1,…,n\)
Varianza constante: \(Var(e_i)=\sigma^2 I\) \(i=1,…,n\)
Residuos incorrelacionados: \(Cov(e_i,e_j) = 0\)
El incumplimiento de alguna de dichas hipótesis, implica la no aleatoriedad de los residuos y, por tanto, la existencia de alguna estructura o relación de dependencia en los residuos que puede ser estimada, debiendo ser considerada en la especificación inicial del modelo. Los principales problemas asociados al incumplimiento de las hipótesis de normalidad de los residuos son, por un lado, la heteroscedasticidad, cuando la varianza de los mismos no es constante, y la autocorrelación o existencia de relación de dependencia o correlación entre los diferentes residuos, lo que violaría el supuesto de términos de error incorrelacionados.
Si se construye una gráfica de los resultados de una estimación mínimo cuadrática (en ordenadas) frente al valor absoluto de los residuos (en abscisas), cuando éstos últimos presentan una distribución aleatoria, es decir una distribución Normal de media cero y varianza constante, , el resultado obtenido (véase Figura 3.1.) muestra que el tamaño del error es independiente del tamaño de la variable estimada, ya que errores con valor elevado se corresponden con valores bajos y altos de la variable dependiente estimada; sin embargo, una distribución de residuos con problemas de heteroscedasticidad da lugar a una figura como la que puede observarse en la Figura 3.2., en donde se manifiesta una clara relación de dependencia entre la variable estimada y el tamaño del error. En este caso los errores de mayor tamaño se corresponden con los valores más altos de la variable estimada.
Figura 3.1 Residuos homocedásticos
Figura 3.2 Residuos heterocedásticos
La representación gráfica de los errores en forma de serie temporal, es decir, poniendo en el eje de ordenadas los errores y en abscisas el periodo temporal en que están datados, permite apreciar la ausencia o presencia de correlación ya que a los residuos no correlacionados (Figura 3.3.) les corresponde una representación gráfica en la que no se aprecia pauta temporal alguna, sucediéndose de forma impredecible o aleatoria, mientras que en los residuos con problemas de autocorrelación la pauta temporal es evidente, evidenciándose que cada residuo podría ser previsto en función de la sucesión de los errores correspondientes a periodos temporales pasados (Figura 3.4.)
Figura 3.3. Residuos sin Autocorrelación
Figura 3.4. Residuos con Autocorrelación
Estos problemas asociados a los errores pueden detectarse con tests estadísticos diseñados para ello.
Heterocedasticidad
La detección de la heteroscedasticidad se realiza a través de diversos contrastes paramétricos, entre los que cabe destacar el contraste de Bartlett (Mood, 1950), el constraste de Goldfeld-Quandt (1965) y el contraste de White (1980).
El contraste de White se basa en que, bajo la hipótesis nula de homocedasticidad, la matriz de varianzas y covarianzas de los estimadores MCO de \(\hat \beta\) es:\(\sigma^2(X'X) ^{-1}\)
Por el contrario, si existe heteroscedasticidad, la matriz de varianzas y covarianzas viene dada por:
\[(X'X) ^{-1}X'\Omega X(X'X)^{-1},\Omega=(\sigma_1^2,\sigma_2^2,...\sigma_n^2) \]
Por tanto, si tomamos la diferencia entre ambas queda:
\[(X'X) ^{-1}X'\Omega X(X'X)^{-1}-\sigma^2(X'X) ^{-1} \]
Por ello, basta con contrastar la hipótesis nula de que todas estas diferencias son iguales a cero, lo que equivale a contrastar que no hay heteroscedasticidad.
Los pasos a seguir para realizar el contraste de White son los siguientes:
Estimar el modelo original y obtener la serie de residuos estimados
Realizar una regresión del cuadrado de la serie de residuos obtenidos en el paso anterior sobre una constante, las variables exógenas del modelo original, sus cuadrados y los productos cruzados de segundo orden (los productos resultantes de multiplicar cada variable exógena por cada una de las restantes). Es decir, se trata de estimar por MCO la relación:
\[e_i^2= \alpha + \varphi_2 X_{2i}+...+\varphi_k X_{ki}+ \eta_2 X_{2i}^2+...+\eta_k X_{ki}^2+ \varpi_3 X_{2i}X_{3i}+...+ \varpi_k X_{2i}X _{ki}+ ....+\rho X_{k-1i }X_{ki} + u_i\]
- Al aumentar el tamaño muestral, el producto \(nR^2\) (donde n es el número de observaciones y \(R^2\) es el coeficiente de determinación de la última regresión) sigue una distribución Chi-cuadrado con \(p – 1\) grados de libertad, donde \(p\) es el número de variables exógenas utilizadas en la segunda regresión. Se aceptará la hipótesis de existencia de heteroscedasticidad cuando el valor del estadístico supere el valor crítico de la distribución Chi-cuadrado (\(c\)) al nivel de significación estadística fijado (\(nR^2>c\)).
Para realizar en R el constraste de heterocedasticidad de White en el modelo estimado para las distancias y velocidades de la base de datos “cars” sobre distancias y velocidades de frenado de automóviles de los años 20 que incorpora R, primero hay que instalar (función “install”) y llamar al Package-R: “tseries”:
library(tseries)
y <- matrix(cars$dist,ncol=1)
x <- matrix(cars$speed,ncol=1)
white.test(x,y)##
## White Neural Network Test
##
## data: x and y
## X-squared = 2.5789, df = 2, p-value = 0.2754Autocorrelación
Decimos que existe autocorrelación cuando el término de error de un modelo econométrico está correlacionado consigo mismo a través del tiempo, tal que \(Cov(e_i,e_j) \neq 0\). Ello no significa que la correlación entre los errores se dé en todos los periodos sino que puede darse tan sólo entre algunos de ellos. En presencia de autocorrelación, los estimadores MCO siguen siendo insesgados pero no poseen mínima varianza, debiéndose utilizar en su lugar el método de estimación de los Mínimos Cuadrados Generalizados (MCG).
La existencia de autocorrelación en los residuos es fácilmente identificable obteniendo las funciones de autocorrelación (acf) y autocorrelación parcial (acp) de los errores mínimo-cuadráticos obtenidos en la estimación. Si dichas funciones corresponden a un ruido blanco, se constatará la ausencia de correlación entre los residuos. Sin embargo, el mero examen visual de las funciones anteriores puede resultar confuso y poco objetivo, por lo que en la práctica econométrica se utilizan diversos contrastes para la autocorrelación, siendo el más utilizado el de Durbin-Watson (1950), que pasamos a ver seguidamente.
Si se sospecha que el término de error del modelo econométrico tiene una estructura como la siguiente:
\[\hat e_t=\rho\hat e_{t-1}+u_t\]
entonces el contraste de Durbin-Watson permite contrastar la hipótesis nula de ausencia de autocorrelación.
Dicho contraste se basa en el cálculo del estadístico \(d\), utilizando para ello los errores mínimo-cuadráticos resultantes de la estimación:
\[d=\frac {\sum_{t=2}^n (\hat e_t - \hat e_{t-1})^2}{\sum_{t=1}^n \hat e_t^2}\]
El valor del estadístico d oscila entre 0 y 4, siendo los valores cercanos a 2 los indicativos de ausencia de autocorrelación de primer orden. La interpretación exacta del test resulta compleja, ya que los valores críticos apropiados para contrastar la hipótesis nula de no autocorrelación requieren del conocimiento de la distribución de probabilidad bajo el supuesto de cumplimiento de dicha hipótesis nula, y dicha distribución depende a su vez de los valores de las variables explicativas, por lo que habría que calcularla en cada aplicación. Se ha demostrado que el valor esperado de d, cuando \(\rho=0\), está dado por la siguiente relación (Maddala,1996).
\[E(d) \approx 2+ \frac {2(k-1)}{n-k}\]
Para facilitar la interpretación del test, Durbin y Watson derivaron dos distribuciones: \(d_L\) y \(d_U\), que no dependen de las variables explicativas y entre las cuales se encuentra la verdadera distribución de \(d\), de forma que a partir de un determinado nivel de significación, se adoptan la siguientes reglas de decisión:
Si \(d \leq d_L\) rechazamos la hipótesis nula de no autocorrelación frente a la hipótesis alternativa de autocorrelación positiva.
Si \(d \geq (4-d_L)\) rechazamos la hipótesis nula de no autocorrelación frente a la hipótesis alternativa de autocorrelación negativa.
Si \(d_U \leq d \leq (4-d_U)\) aceptamos la hipótesis nula de no autocorrelación.
Si \(d_L \leq d \leq d_U\) ó \((4-d_U) \leq d \leq (4-d_L)\) la prueba no es concluyente.
La Figura siguiente muestra esta interpretación de forma gráfica:
Interpretación test de Durbin-Watson
El estadístico de Durbin-Watson es aproximadamente igual a \(2(1-r)\) en donde \(r\) es el coeficiente de autocorrelación simple muestral del retardo 1.
En R, el test de Durbin-Watson se encuentra en el Package-R: “lmtest”, y su sintaxis es:
library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericdwtest(y~x)##
## Durbin-Watson test
##
## data: y ~ x
## DW = 1.6762, p-value = 0.09522
## alternative hypothesis: true autocorrelation is greater than 0Deficiencias muestrales
El fenómeno de la multicolinealidad aparece cuando las variables exógenas de un modelo econométrico están correlacionadas entre sí, lo que tiene consecuencias negativas para la estimación por MCO, ya que la existencia de una relación lineal entre las variables exógenas, implica que la matriz (\(X'X\)) va a tener determinante cero, es decir será una matriz singular y por tanto no será invertible. Dado que \(\hat \beta = X'X ^{-1} X'y\), no será posible calcular la estimación mínimo cuadrática de los parámetros del modelo ni, lógicamente, la varianza de los mismos. Esto es lo que se conoce por el nombre de multicolinealidad exacta.
Errores de especificación
Los errores de especificación hacen referencia a un conjunto de errores asociados a la especificación de un modelo econométrico. En concreto cabe referirse a:
- Omisión de variables relevantes
- Inclusión de variables innecesarias
- Adopción de formas funcionales equivocadas
En Economía la teoría no suele concretar la forma funcional de las relaciones que estudia. Así, por ejemplo, cuando se analiza la demanda se señala que la cantidad demandada es inversamente proporcional al precio; cuando se estudia el consumo agregado se apunta que la propensión marginal a consumir (relación entre renta y/o consumo) es mayor que cero y menor que uno. Por otro lado es frecuente utilizar la condición “ceteris paribus” para aislar la información de otras variables relevantes que influyen y/o modifican la relación estudiada. Por esta razón, la existencia de errores de especificación en la relación estimada es un factor a considerar y a resolver en el proceso de la estimación econométrica.
Con independencia de la naturaleza de los errores de especificación, dado que el proceso de estimación MCO deben de cumplirse determinadas hipótesis básicas, que los estimadores MCO deben de ser insesgados, eficientes y consistentes, y que el estimador de la varianza del término de error ha de ser insesgado, debemos preguntarnos: ¿qué ocurriría con estas propiedades ante errores de especificación?.
La omisión de variables relevantes provoca que el estimador MCO sea sesgado en media (se obtiene un valor distinto que el que se obtendría al incorporar la variable relevante omitida) y en varianza. La incorporación de variables innecesarias determina que la varianza de los estimadores MCO sea mas elevada (sesgo en la varianza), dificultando la interpretación del contraste de significación individual sobre los parámetros.
Si especificamos la forma funcional de una relación (ya sea lineal, cuadrática, cúbica, exponencial, logarítmica, etc.) y la verdadera relación presenta una forma diferente a la especificada tiene, en algunos casos, las mismas consecuencias que la omisión de variables relevantes, es decir, proporciona estimadores sesgados e inconsistentes. En general, una especificación funcional incorrecta lleva a obtener perturbaciones heteroscedásticas y/o autocorrelacionadas, o alejadas de los parámetros de la distribución del término de error del modelo correctamente especificado.
Métodos de selección de variables en el modelo lineal general
Una de las cuestiones más importantes a la hora de encontrar el modelo de ajuste más adecuado cuando se dispone de un amplio conjunto de variables explicativas, es la correcta especificación del modelo teórico, ya que como se ha visto la inclusión de una variable innecesaria o la omisión de una variable relevante, condiciona los estadísticos que resultan en la estimación MCO del modelo, por otro lado, en un elevado número de explicativas no cabe descartar la existencia de correlaciones que originen un problema de multicolinealidad aproximada, y en estos casos hay que determinar cual de ellas cabe incluir en la especificación del modelo.
En otras palabras, ante un conjunto elevado de explicativas debemos seleccionar de entre todas, un subconjunto de ellas que garanticen que el modelo esté lo mejor especificado posible. Este análisis cabe hacerlo estudiando las características y propiedades de cada una de las variables independientes, a partir, por ejemplo, de los coeficientes de correlación de cada una de ellas y la dependiente, y de cada explicativa con las restantes, seleccionando modelos alternativos y observando los resultados estadísticos de la estimación MCO de cada uno de ellos. Sin embargo, en la práctica, la selección del subconjunto de variables explicativas de los modelos de regresión se deja en manos de procedimientos más o menos automáticos.
Los procedimientos más usuales son los siguientes:
• Método backward: se comienza por considerar incluidas en el modelo teórico a todas las variables disponibles y se van eliminando del modelo de una en una según su capacidad explicativa. En concreto, la primera variable que se elimina es aquella que presenta un menor coeficiente de correlación parcial con la variable dependiente-o lo que es equivalente, un menor valor del estadístico t– y así sucesivamente hasta llegar a una situación en la que la eliminación de una variable más suponga un descenso demasiado acusado en el coeficiente de determinación.
• Método forward: se comienza por un modelo que no contiene ninguna variable explicativa y se añade como primera de ellas a la que presente un mayor coeficiente de correlación -en valor absoluto- con la variable dependiente. En los pasos sucesivos se va incorporando al modelo aquella variable que presenta un mayor coeficiente de correlación parcial con la variable dependiente dadas las independientes ya incluidas en el modelo. El procedimiento se detiene cuando el incremento en el coeficiente de determinación debido a la inclusión de una nueva variable explicativa en el modelo ya no es importante.
• Método stepwise: es uno de los más empleados y consiste en una combinación de los dos anteriores. En el primer paso se procede como en el método forward pero a diferencia de éste, en el que cuando una variable entra en el modelo ya no vuelve a salir, en el procedimiento stepwise es posible que la inclusión de una nueva variable haga que otra que ya estaba en el modelo resulte redundante.
El modelo de ajuste al que se llega partiendo del mismo conjunto de variables explicativas es distinto según cuál sea el método de selección de variables elegido, por lo que la utilización de un procedimiento automático de selección de variables no significa que con él se llegue a obtener el mejor de los modelos a que da lugar el conjunto de datos con el que se trabaja.
Para realizar la selección de un modelo por cualquira de los métodos descritos, necesitamos instalar la librería-R: “leaps”, una vez instalada ejecutamos el siguiente Chunk, para seleccionar según el método “forward” :
library(leaps)
regfit.fwd=regsubsets(mpg~.,data=mtcars,method="forward")
plot(regfit.fwd,scale="r2")
summary(regfit.fwd) ## Subset selection object
## Call: regsubsets.formula(mpg ~ ., data = mtcars, method = "forward")
## 10 Variables (and intercept)
## Forced in Forced out
## cyl FALSE FALSE
## disp FALSE FALSE
## hp FALSE FALSE
## drat FALSE FALSE
## wt FALSE FALSE
## qsec FALSE FALSE
## vs FALSE FALSE
## am FALSE FALSE
## gear FALSE FALSE
## carb FALSE FALSE
## 1 subsets of each size up to 8
## Selection Algorithm: forward
## cyl disp hp drat wt qsec vs am gear carb
## 1 ( 1 ) " " " " " " " " "*" " " " " " " " " " "
## 2 ( 1 ) "*" " " " " " " "*" " " " " " " " " " "
## 3 ( 1 ) "*" " " "*" " " "*" " " " " " " " " " "
## 4 ( 1 ) "*" " " "*" " " "*" " " " " "*" " " " "
## 5 ( 1 ) "*" " " "*" " " "*" "*" " " "*" " " " "
## 6 ( 1 ) "*" "*" "*" " " "*" "*" " " "*" " " " "
## 7 ( 1 ) "*" "*" "*" "*" "*" "*" " " "*" " " " "
## 8 ( 1 ) "*" "*" "*" "*" "*" "*" " " "*" "*" " "Si queremos ver las estimaciones MCO de los parámetros del modelo 8:
coef(regfit.fwd,8)## (Intercept) cyl disp hp drat wt
## 12.56350226 -0.23126963 0.01611609 -0.02339020 0.70893592 -4.08155351
## qsec am gear
## 0.91812006 2.47759723 0.50403957Regresión Ridge
Uno de los supuestos básicos del métodos de mínimos cuadrados es que el rango de \(X\), es decir el número de columnas (o filas) linealmente independientes tiene que ser \(k\), ya que \(\rho(X)=k\), entonces \(\rho(X′X) =k\) y el sistema de ecuaciones normales tiene solución única. Si el supuesto se incumple, \(\rho(X) < k\), entonces las columnas de la matriz \(X\) son linealmente dependientes, \(\rho(X′X) < k\) y el sistema de ecuaciones normales tiene soluciones múltiples. El término multicolinealidad hace referencia a la existencia de una o más relaciones lineales exactas o perfectas entre las variables explicativas.
Este es el caso de la multicolienealidad exacta, si entre los regresores hay variables correlacionadas, estamos ante un problema de multicolinualidad imperfecta, que conduce a sesgos en la estimación de los parámetros asociados a las variables correlacionadas.
La regresión contraída o regresión ridge trata de dar solución a estos problemas de presencia de multicolinealidad, estimando unos parametros “\(\beta\) contraidos”:
\(\hat \beta^k = (X’X + \lambda I)^{-1}X’y\)
Si lambda es 0 estamos ante mínimos cuadrados ordinarios, en otro caso estamos ante un estimador sesgado de \(\beta\). Este estimador sesgado es la solución a un problema de mínimos cuadrados penalizados que contrae los \(\beta\) en torno a 0.
En definitiva, se trata de de encontrar el coeficiente \(\lambda\) que contrae las estimaciones de \(\beta\), y para encontrarlo se utiliza un criterio generalizado de validación cruzada, generalized cross-validation (GCV). El proceso fija un rango de posibles \(\lambda\) entre \([0, c]\) y calcula la validación cruzada \(CV\), siendo el \(\lambda\) óptimo aquel que minimiza el CV.
A continuación se presenta un ejemplo de un modelo \(Y_{t} = \beta_0 + \beta_1 X_{1,t} + \beta_2 X_{2,t} + + \beta_3 X{3,t} e_t\). Donde las variables son generadas de forma aleatoria, en \(X_{3c}\) se introduce una dependencia lineal entre las variables \(X_1\) y \(X_3\).
La estimación de la regresión ridge se realiza con la función “lm.ridge” de la liberia “MASS”, fijando un rango de posibles \(\lambda\) entre \([0, 10]\). Se elige finalmente un \(\lambda=4\).
# Three variables are measured: x1,x2,x3.
# x1 and x1 are U(0,1); x3=10 * X1 + unif(0,1).
# This causes corr(X1,X3)=sqrt(100/101)=0.995.
# We will fit OLS and ridge regressions to these data,
# use the data to select the "best" constant to add,
# and then evaluate the two regressions on a new test set.
# Ridge regression function, ridge.lm(), is on MASS package
library(MASS)
# Generating the data
set.seed(558562316)
N <- 20 # Sample size
x1 <- runif(n=N)
x2 <- runif(n=N)
x3 <- runif(n=N)
x3c <- 10*x1 + x3 # New variable
ep <- rnorm(n=N)
y <- x1 + x2 + ep
# OLS fit of 3-variable model using independent x3
ols <- lm(y~ x1 + x2 + x3)
summary(ols)##
## Call:
## lm(formula = y ~ x1 + x2 + x3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.19698 -0.28592 0.04026 0.24016 1.20322
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.4293 0.4916 -0.873 0.3954
## x1 1.7851 0.4812 3.710 0.0019 **
## x2 0.7119 0.4622 1.540 0.1430
## x3 0.2839 0.5122 0.554 0.5870
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6306 on 16 degrees of freedom
## Multiple R-squared: 0.4831, Adjusted R-squared: 0.3862
## F-statistic: 4.984 on 3 and 16 DF, p-value: 0.0125# OLS fit of 3-variable model using correlated x3.
olsc <- lm(y~ x1 + x2 + x3c)
summary(olsc)##
## Call:
## lm(formula = y ~ x1 + x2 + x3c)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.19698 -0.28592 0.04026 0.24016 1.20322
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.4293 0.4916 -0.873 0.395
## x1 -1.0540 5.1966 -0.203 0.842
## x2 0.7119 0.4622 1.540 0.143
## x3c 0.2839 0.5122 0.554 0.587
##
## Residual standard error: 0.6306 on 16 degrees of freedom
## Multiple R-squared: 0.4831, Adjusted R-squared: 0.3862
## F-statistic: 4.984 on 3 and 16 DF, p-value: 0.0125# Ridge regression using independent variables
ridge <- lm.ridge (y ~ x1+x2+x3, lambda = seq(0, .1, .001))
summary(ridge)## Length Class Mode
## coef 303 -none- numeric
## scales 3 -none- numeric
## Inter 1 -none- numeric
## lambda 101 -none- numeric
## ym 1 -none- numeric
## xm 3 -none- numeric
## GCV 101 -none- numeric
## kHKB 1 -none- numeric
## kLW 1 -none- numericplot(ridge)
# Ridge regression using correlated variables
ridgec <- lm.ridge (y ~ x1+x2+x3c, lambda = seq(0, .1, .001))
plot(ridgec)
select(ridgec)## modified HKB estimator is 0.4209458
## modified L-W estimator is 1.337588
## smallest value of GCV at 0.1# Selection of constant is at endpoint. Extend endpoint and try again
ridgec <- lm.ridge (y ~ x1+x2+x3c, lambda = seq(0, 1, .1))
plot(ridgec)
select(ridgec)## modified HKB estimator is 0.4209458
## modified L-W estimator is 1.337588
## smallest value of GCV at 1# Selection of constant is at endpoint. Extend endpoint and try again
ridgec <- lm.ridge (y ~ x1+x2+x3c, lambda = seq(0, 10, .01))
plot(ridgec)
select(ridgec)## modified HKB estimator is 0.4209458
## modified L-W estimator is 1.337588
## smallest value of GCV at 4.01# Final model uses lambda=4.
ridge.final <- lm.ridge (y ~ x1+x2+x3c, lambda = 4)
ridge.final## x1 x2 x3c
## -0.13809086 0.76302466 0.56030303 0.08372054summary(ridge.final)## Length Class Mode
## coef 3 -none- numeric
## scales 3 -none- numeric
## Inter 1 -none- numeric
## lambda 1 -none- numeric
## ym 1 -none- numeric
## xm 3 -none- numeric
## GCV 1 -none- numeric
## kHKB 1 -none- numeric
## kLW 1 -none- numeric# Create test data and compute predicted values for OLS and ridge.
# There's no predict() method for "ridgelm" objects
test <- expand.grid(x1 = seq(.05,.95,.1), x2 = seq(.05,.95,.1), x3=seq(.05,.95,.1))
mu <- test$x1 + test$x2
test$x3c <- 10*test$x1 + test$x3
pred.ols <- predict(ols,newdata=test) # y ~ X1 + X2 + X3
pred.olsc <- predict(olsc,newdata=test) # y ~ X1 + X2 + X3c
pred.ridge <- coef(ridge.final)[1] + coef(ridge.final)[2]*test[,1] +
coef(ridge.final)[3]*test[,2] + coef(ridge.final)[4]*test[,4]
MSPE.ols <- sum((pred.ols - mu)^2)/length(mu)
MSPE.olsc <- sum((pred.olsc - mu)^2)/length(mu)
MSPE.ridge <- sum((pred.ridge - mu)^2)/length(mu)
MSPE.ols## [1] 0.06585964MSPE.olsc## [1] 0.06585964MSPE.ridge## [1] 0.04650589