Estadística y Machine Learning con R

Portada

Introducción

El manual de Estadística y Machine Learning con R es el resultado de diversas experiencias que ha tenido lugar en el Instituto Cantabro de Estadística (ICANE), en lo relativo a la formación del personal. Una vez se considero utilizar el software estadístico R, como referente para los futuros desarrollos estadísticos, sutituyendo a diversos software (SPSS, Gretl, Octave, etc…) que utilizaban los diferentes departamentos, se desarrolló un plan de formación en R, que se inició en Octubre de 2016, con el curso “Introducción al análisis de datos con R”, impartido por los profesores Miguel Expósito Martín y Francisco Parra Rodríguez, dentro de la programación de cursos para formación para el gobierno de Cantabria del CEAR. Ese mismo año se impartió un curso de tecnicas estadísticas en R para el personal funcionario y becario del ICANE. Para ámbos cursos se elaboró un documento de trabajo (DOC nº2 2016): Curso de Estadística con R, en donde junto a las técnicas estadísticas se introducían algunos elementos del aprendizaje automático (http://www.icane.es/c/document_library/get_file?uuid=c2e9fff1-72d5-42ad-b391-bacb3ebe9dbe&groupId=10138). En 2017, se impartió a los becarios otro curso de formación en R, en donde se ampliaron los contenidos de machine learning y los de estadística multivariante. El propósito de dicha formación era que las personas con beca en el ICANE, dispusieran habilidades e instrumentos de libre difusión con los que pudieran completar su capacitación estadística en el futuro. Dado que el perfil de las becas es amplio en lo que a titulaciones que permite (estadísticos, economístas, demógrafos, sociologos, etc…) se planteó una formación de caracter lo más generalista posible y no centrada en la producción estadística que realiza el ICANE.

El manual de Estadística y Machine Learnig con R, a diferencia de otros documentos técnicos del ICANE, ha sido elaborado exclusivamente con R, utilizando las librerias knitr, markdown y bookdown, que permiten editar y compilar documentos en diferentes formatos. Con el fin de que pueda utilizarse en futuros proyectos formativos, el documento se ha compilado en html, para ser consultado on-line y difundido a través de las redes sociales.

Modelo Lineal General

Introducción

La regresión lineal es la técnica básica del análisis econométrico. Mediante dicha técnica tratamos de determinar relaciones de dependencia de tipo lineal entre una variable dependiente o endógena, respecto de una o varias variables explicativas o exógenas. Gujarati (1975), define el análisis de regresión como el estudio de la dependencia de la variable dependiente, sobre una o más variables explicativas, con el objeto de estimar o predecir el valor promedio poblacional de la primera en términos de los valores conocidos o fijos (en medias muestrales repetidas) de las últimas.

En este capitulo abordaremos el estudio del caso de una única ecuación de tipo lineal con una variable dependiente y una independiente, y la generalización del modelo al caso de múltiples variables exógenas. Las extensiones del modelo lineal general se analizarán en capítulos siguientes.

Modelo de Regresión Lineal Simple

Partimos de la existencia de una relación lineal entre una variable endógena (\(Y_i\)) y una variable exógenas (\(X_i\)):

\[\begin{equation} Y_i= \beta_1 + \beta_2 X_i + e_i, i=1,2,…, n \end{equation}\]

Nuestro objetivo consiste en estimar los dos parámetros de la ecuación anterior a partir de los datos muestrales de los que disponemos. Para ello utilizaremos el método de los Mínimos Cuadrados Ordinarios (MCO), pero antes de ver en que consiste este método debemos plantear ciertas hipótesis sobre el comportamiento de las variables que integran el modelo.

La variable \(e_i\) la denominamos término de perturbación o error, y en ella recogemos todos aquellos factores que pueden influir a la hora de explicar el comportamiento de la variable \(Y_i\) y que, sin embargo, no están reflejados en las variables explicativas, \(X_i\). Estos factores deberían ser poco importantes, ya que no debería existir ninguna variable explicativa relevante omitida en el modelo de regresión. En caso contrario estaríamos incurriendo en lo que se conoce como un error de especificación del modelo. El término de perturbación también recogería los posibles errores de medida de la variable dependiente, \(Y_i\).

De lo anterior se desprende que, a la hora de estimar los parámetros del modelo, resultará de vital importancia que dicho término de error no ejerza ninguna influencia determinante en la explicación del comportamiento de la variable dependiente. Por ello, si el modelo está bien especificado, cuando se aplica el método de Mínimos Cuadrados Ordinarios, cabe realizar las siguientes hipótesis de comportamiento sobre el término de error:

1. La esperanza matemática de et es cero, tal que \(E(e_i) = 0\). Es decir, el comportamiento del término de error no presenta un sesgo sistemático en ninguna dirección determinada. Por ejemplo, si estamos realizando un experimento en el cual tenemos que medir la longitud de un determinado objeto, a veces al medir dicha longitud cometeremos un error de medida por exceso y otras por defecto, pero en media los errores estarán compensados.

2. La covarianza entre \(e_i\) y \(e_j\) es nula para tal que \(E (e_i•e_j) = 0\). Ello quiere decir que el error cometido en un momento determinado, i, no debe estar correlacionado con el error cometido en otro momento del tiempo, j, o dicho de otro modo, los errores no ejercen influencia unos sobre otros. En caso de existir este tipo de influencia o correlación, nos encontraríamos ante el problema de la autocorrelación en los residuos, el cual impide realizar una estimación por Mínimos Cuadrados válida.

3. La matriz de varianzas y covarianzas del término de error debe ser escalar tal que \(Var(e_i)=\sigma^2I\), \(i=1,…,n\), donde I es la matriz unidad. Dado que siempre que medimos una variable, se produce un cierto error, resulta deseable que los errores que cometamos en momentos diferentes del tiempo sean similares en cuantía. Esta condición es lo que se conoce como supuesto de homocedasticidad que, en caso de no verificarse, impediría un uso correcto de la estimación lineal por Mínimos Cuadrados.

Estas hipótesis implican que los errores siguen una distribución Normal de media cero y varianza constante por lo que, dado su carácter aleatorio, hace que los errores sean por naturaleza impredecibles.

Asimismo, las variables incluidas en el modelo deben verificar que:

1. La variable dependiente \(Y_i\) se ajusta al modelo lineal durante todo el periodo muestral, es decir, no se produce un cambio importante en la estructura de comportamiento de \(Y_i\) a lo largo de la muestra considerada, es decir, se distribuira como una normal con \(\mathrm{E}((Y_i)=\beta_1 + \beta_2 X_i\) y \(\mathrm{Var}(Y_i)=\mathrm{Var}(e_i)\).

2. La variable explicativa, \(X_i\), es no estocásticas, es decir, es considerada fija en muestreos repetidos.

Si suponemos que se verifican los supuestos anteriores, la estimación mínimo cuadrática de los parámetros \(\beta_1\) y \(\beta_2\), dará como resultado gráfico una recta que se ajuste lo máximo posible a la nube de puntos definida por todos los pares de valores muestrales \((X_i,Y_i)\), tal y como se puede apreciar en el Figura 2.1.

Figura 2.1

El término de error, \(e_i\), puede ser entendido, a la vista del gráfico anterior, como la distancia que existe entre el valor observado, \(Y_i\), y el correspondiente valor estimado, que sería la imagen de \(X_i\) en el eje de ordenadas. El objetivo de la estimación por Mínimos Cuadrados Ordinarios es, precisamente, minimizar el sumatorio de todas esas distancias al cuadrado; es decir :

\[\begin{equation} Min \sum (e_i)^2= \sum (Y_i - \hat Y_i)^2= \sum (Y_i - \hat\beta_1 + \hat\beta_2 X_i)^2 \end{equation}\]

Derivando esta expresión respecto a los coeficientes \(\hat\beta_1\) y \(\hat\beta_2\) e igualando a cero obtenemos el sistema de ecuaciones normales:

\[\begin{equation} \sum Y_i= n\hat\beta_1 + \sum \hat\beta_2 X_i \\ \sum (Y_i•X_i)= \sum \hat\beta_1 X_i + \sum \hat\beta_2 (X_i)^2 \end{equation}\]

donde n representa el tamaño muestral.

Resolviendo dicho sistema de ecuaciones obtenemos la solución para los parámetros \(\hat\beta_1\) y \(\hat\beta_2\): \[\begin{equation} \hat\beta_1 =\bar X - \sum \hat\beta_2 \bar X \\ \hat\beta_2=\frac{\sum ((Y_i-\bar Y)•(X_i-\bar X)} {\sum (X_i-\bar X)^2} \end{equation}\]

Modelo de Regresión Lineal Múltiple

Pasamos a continuación a generalizar el modelo anterior al caso de un modelo con varias variables exógenas, de tal forma que se trata de determinar la relación que existe entre la variable endógena Y y variables exógenas: \(X_2, X_3,…, X_k\). Dicho modelo se puede formular matricialmente de la siguiente manera:

\[\begin{equation} y=\beta X+e=Y_i= \beta_1 + \beta_2 X_{2i} +\beta_3 X_{3i} +...+\beta_k X_{ki} e_i, i=1,2,…, n \end{equation}\]

donde: \[y = \begin{bmatrix} Y_{1}\\ Y_{2}\\ .\\ Y_{n} \end{bmatrix}\] es el vector de observaciones de la variable endógena

\[X = \begin{bmatrix}1 & X_{21} & X_{31} & ... & X_{k1}\\ 1 & X_{22} & X_{32} & ... & X_{k2}\\ .& . & . & ... & .\\ 1 & X_{2n} & X_{3n} & ... & X_{kn} \end{bmatrix}\] es el vector de observaciones de las variables exogenas. \[\beta = \begin{bmatrix} \beta_{1}\\ \beta_{2}\\ .\\ \beta_{k} \end{bmatrix}\] es el vector de los coeficientes que pretendemos estimar.

\[e = \begin{bmatrix} e_{1}\\ e_{2}\\ .\\ e_{n} \end{bmatrix}\] es el vector de términos de error.

Suponiendo que se verifican las hipótesis que veíamos antes, el problema a resolver nuevamente es la minimización de la suma de los cuadrados de los términos de error tal que:

\[\begin{equation} Min \sum (e_i)^2= \sum (Y_i - \hat Y_i)^2 =\sum (Y_i - \hat\beta X)^2 \end{equation}\]

Desarrollando dicho cuadrado y derivando respecto a cada \(\hat\beta_i\) obtenemos el siguiente sistema de ecuaciones normales expresado en notación matricial:

\[ X'X\hat \beta = X'y\]

en donde basta con despejar \(\hat\beta\) premultiplicando ambos miembros por la inversa de la matriz para obtener la estimación de los parámetros del modelo tal que:

\[ \hat \beta =X'X ^{-1} X'y\]

donde:

\[X'X = \begin{bmatrix}n & \sum X_{2i} & \sum X_{3i} & ... & \sum X_{ki}\\ \sum X_{2i} & \sum X_{2i}^2 & \sum (X_{2i}X_{3i}) & ... & \sum (X_{2i}X_{ki})\\ .& . & . & ... & .\\ \sum X_{ki} & \sum (X_{ki}X_{2i}) & \sum (X_{ki}X_{3i}) & ... & \sum X_{ki}^2 \end{bmatrix}\]

\[X'y = \begin{bmatrix} \sum Y_{i}\\ \sum (Y_{i}X_{2i})\\ \\ \sum (Y_{i}X_{ki}) \end{bmatrix}\]

\[\hat \beta = \begin{bmatrix} \hat \beta_{1}\\ \hat \beta_{2}\\ .\\ \hat \beta_{k} \end{bmatrix}\]

Cada uno de los coeficientes estimados,\(\hat \beta_i\), son una estimación insesgada del verdadero parámetro del modelo y representa la variación que experimenta la variable dependiente \(Y_i\) cuando una variable independiente \(X_i\) varía en una unidad y todas las demás permanecen constantes (supuesto ceteris paribus). Dichos coeficientes poseen propiedades estadísticas muy interesantes ya que, si se verifican los supuestos antes comentados, son insesgados, eficientes y óptimos.

La valores estimados de la variable dependiente en expresión matricial vienen dados por:

\[\hat y = X \hat \beta \]

Y la distribución estadística de la variable dependiente, se expresa en forma matricial:

\[\mathrm{E}(y)= X \beta\]

\[\mathrm{Var}(y)=\sigma^{2}I\]

Propiedades estadísticas del estimador MCO

Si se cumplen las hipótesis de comportamiento sobre el término error, la distribución de probabilidad del estimador MCO \(\hat \beta\) será una distribución Normal multivariante con vector de medias \(\beta\) y matriz de varianzas y covarianzas \(\sigma^2(X'X) ^{-1}\). La esperanza matemática del estimador MCO se demuestra a partir de:

\[ \mathrm{E}(\hat \beta) = \mathrm{E}(X'X ^{-1}X'u) \\ =X'X ^{-1}X'\mathrm{E}(y) \\ =X'X^{-1}X'X\beta= \beta \]

De la definición de matriz de varianzas y covarianzas, y se tiene que:

\[\begin{equation} \begin{split} \mathrm{Var}(\hat{\beta}) & =\mathrm{Var}((X'X)^{-1}X'y)\\ & =(X'X)^{-1}X'\mathrm{Var}(y)((X'X)^{-1}X')'\\ & =(X'X)^{-1}X'\mathrm{Var}(y)X(X'X)^{-1}\\ & =(X'X)^{-1}X'\sigma^{2}IX(X'X)^{-1}\\ & =(X'X)^{-1}\sigma^{2} \end{split} \end{equation}\]

El estimador \(\hat \beta\) del parámetro \(\beta\) es insesgado porque su esperanza matemática coincide con el verdadero valor del parámetro \(\mathrm{Var}(\hat \beta)=\beta\) . Se dice que un estimador insesgado es mas eficiente que otro estimador insesgado, si la varianza muestral de dicho estimado res menor que la varianza muestral de ese otro estimador. El teorema de Gauss-Markov demuestra que el estimador MCO es el más eficiente de la clase de estimadores lineales e insesgados de \(\beta\).

Coeficiente de determinación

Una vez estimada la ecuación de regresión lineal tiene interés determinar la exactitud del ajuste realizado. Para ello hay que analizar la variación que experimenta esta variable dependiente y, dentro de esta variación, se estudia qué parte está siendo explicada por el modelo de regresión y qué parte es debida a los errores o residuos.

La forma de realizar dicho análisis es a partir de la siguiente expresión:

\[SCT=SCE+SCR\]

donde:

• \(SCT\): es la Suma de Cuadrados Totales y representa una medida de la variación de la variable dependiente.

• \(SCE\) es la Suma de Cuadrados Explicados por el modelo de regresión.

• \(SCR\) es la Suma de Cuadrados de los Errores

Cuando el modelo tiene término independiente, cada una de estas sumas viene dada por:

\[SCT=y'y-n\bar Y^2=\sum Y_i^2 -n\bar Y^2\]

\[SCT=\hat \beta' X'y \hat \beta-n\bar Y^2=\sum \hat Y_i^2 -n\bar Y^2\]

\[SCE=y'y - \hat \beta' X'y=\sum \hat e_i^2\]

A partir de las expresiones anteriores es posible obtener una medida estadística acerca de la bondad de ajuste del modelo mediante lo que se conoce como coeficiente de determinación (\(R^2\)), que se define como:

\[R^2=\frac{SCE}{SCT}=1-\frac{SCR}{SCT}, 0 \leq R^2 \leq 1\]

Mediante este coeficiente es posible seleccionar el mejor modelo de entre varios que tengan el mismo número de variables exógenas, ya que la capacidad explicativa de un modelo es mayor cuanto más elevado sea el valor que tome este coeficiente. Sin embargo, hay que tener cierto cuidado a la hora de trabajar con modelos que presenten un R2 muy cercano a 1 pues, aunque podría parecer que estamos ante el modelo “perfecto”, en realidad podría encubrir ciertos problemas de índole estadística como la multicolinealidad que veremos en el capítulo 3.

Por otra parte, el valor del coeficiente de determinación aumenta con el número de variables exógenas del modelo por lo que, si los modelos que se comparan tienen distinto número de variables exógenas, no puede establecerse comparación entre sus \(R^2\). En este caso debe emplearse el coeficiente de determinación corregido (\(\bar R^2\)) , el cual depura el incremento que experimenta el coeficiente de determinación cuando el número de variables exógenas es mayor.

La expresión analítica de la versión corregida es:

\[R^2=1-\frac{\frac{SCR}{n-k}}{\frac{SCT}{n-1}}\]

Inferencia acerca de los estimadores

Hasta el momento hemos visto como la estimación por MCO permite obtener estimaciones puntuales de los parámetros del modelo. La inferencia acerca de los mismos permite completar dicha estimación puntual, mediante la estimación por intervalos y los contrastes de hipótesis. Los primeros posibilitan la obtención de un intervalo dentro del cual, con un determinado nivel de confianza, oscilará el verdadero valor de un parámetro, mientras que los segundos nos permitirán extraer consecuencias del modelo, averiguando si existe o no, evidencia acerca de una serie de conjeturas que pueden plantearse sobre sus parámetros.

La inferencia estadística consiste en la estimación de los parámetros poblacionales a partir de la información extraída de una muestra de dicha población. El número de estimaciones que podemos realizar de una población, a través de la extracción de diferentes muestras de un mismo tamaño, es generalmente muy grande porque cada una de las muestras posibles que se pueden sacar de la población arrojaría una estimación.

Por esta razón, a la estimación que obtenemos en una investigación por muestreo la acompañamos con un intervalo de valores posibles. La amplitud de dicho intervalo dependerá del grado de confianza que establezcamos.

El grado o nivel de confianza nos expresa el número de veces que la media verdadera de la población está incluida en cien intervalos de cien muestras extraídas de una población dada. El nivel de confianza más utilizado es el 95%, lo que quiere decir que 95 de cada 100 intervalos construidos contendrán el verdadero valor de la media.

El intervalo de confianza para la media de una población normalmente distribuida se construye en base a la probabilidad de que dicha media esté comprendida entre dos valores \(\bar X_a\) y \(\bar X_b\) equidistantes a ella:

\[\mathrm{P}(\bar X_a \leq \mu_{\bar X} \leq \bar X_b)=1-\alpha\]

siendo \(1-\alpha\) el nivel o grado de confianza asociado a dicho intervalo.

En términos generales, los intervalos de confianza para los estadísticos muestrales se expresan como:

Estimador ± (Factor de Fiabilidad)(Error Típico del Estimador)

Intervalos de Confianza

Para construir los intervalos de confianza de los parámetros \(\beta_i\), se parte de que la estimación MCO proporciona el valor medio de los posibles valores que pudiera tener dicho parámetro, y de que la distribución de dichos valores sigue una distribución derivada de la Normal que se conoce como t de Student.

Dicha distribución es simétrica presentando mayor dispersión que la curva Normal estándar para un tamaño muestral \(n\) pequeño. A medida que \(n\) aumenta (\(n > 100\)) es prácticamente igual que la distribución Normal.

El cálculo del intervalo de confianza para se realiza mediante la siguiente expresión:

\[\hat \beta_i \pm S_{\hat \beta_i}t_{n-k}\]

donde \(S_{\hat \beta_i}\) es la desviación típica estimada para el coeficiente, que se obtiene de la matriz de varianzas y covarianzas de los estimadores MCO \(\hat Var(\hat \beta) =(X'X)^{-1} \hat \sigma^{2}\), calculado \(\hat \sigma^{2}\) a partir de:

\[\hat \sigma^{2}=\frac{\sum \hat e_i^2}{n-k}=\frac{e'e}{n-k}\]

Contraste de hipótesis

Una buena parte de las investigaciones estadísticas están orientadas al desarrollo de procesos encaminados a la contrastación de hipótesis que previamente se han establecido.

Una hipótesis es una afirmación que está sujeta a verificación o comprobación. Hay que tener presente que una hipótesis no es un hecho establecido o firme, las hipótesis están basadas en la experiencia, en la observación, en la experimentación o en la intuición del sujeto que las formula.

Cuando las hipótesis se plantean de tal modo que se pueden comprobar por medio de métodos estadísticos reciben el nombre de hipótesis estadísticas. Estas hipótesis son afirmaciones que se efectúan sobre uno o más parámetros de una o más poblaciones. Las hipótesis estadísticas son de dos tipos: hipótesis nula e hipótesis alternativa. La hipótesis nula, o que no se verifique dicha afirmación, simbolizada por \(H_0\), es la hipótesis que se debe comprobar.

Para contrastar una hipótesis nula examinamos los datos de la muestra tomados de la población y determinamos si son o no compatibles con dicha hipótesis. Si son compatibles entonces \(H_0\) se acepta, en caso contrario se rechaza. Si se acepta la hipótesis nula afirmamos que los datos de esa muestra en concreto no dan suficiente evidencia para que concluyamos que la hipótesis nula sea falsa; si se rechaza decimos que los datos particulares de la muestra ponen de manifiesto que la hipótesis nula es falsa, entonces la hipótesis alternativa,\(H_1\), es considerada verdadera.

El criterio que permite decidir si rechazamos o no la hipótesis nula es siempre el mismo. Definimos un estadístico de prueba, y unos límites que dividen el espacio muestral en una región en donde se rechaza la hipótesis establecida, y otra región en la que no se rechaza, llamada región de aceptación. A la región donde se rechaza la hipótesis nula se le llama región crítica. Esta región es un subconjunto del espacio muestral, y si el valor del estadístico de prueba pertenece a él se rechaza la hipótesis nula.

El límite entre la región crítica y la región de aceptación viene determinado por la información previa relativa a la distribución del estadístico de prueba.

Señalar que un estadístico de prueba es una fórmula que nos dice como confrontar la hipótesis nula con la información de la muestra y es, por tanto, una variable aleatoria cuyo valor cambia de muestra a muestra.

Otra de las consideraciones a realizar en el contraste de hipótesis es fijar la probabilidad del error de rechazar la prueba siendo cierta, a este error se le denomina nivel de significación. Por ejemplo, si se utiliza un nivel de significación de \(\alpha =0,05\), equivale a decir que si para realizar un contraste tomáramos infinitas muestras de la población, rechazaríamos la hipótesis nula de forma incorrecta un 5 % de las veces.

En la formalización del procedimiento de contrastación podemos distinguir siete pasos principales:

1. Formular las hipótesis \(H_0\) y \(H_1\) a contrastar;
2. Proponer un estadístico de contraste;
3. Elegir un nivel de significación aceptable;
4. Determinar la región crítica y/o calcular el p-valor;
5. Tomar la decisión de aceptar o rechazar la hipótesis nula.

Para contrastar la hipótesis de que el parámetro poblacional \(\beta_i\) tome el valor \(\beta_i^0\), necesitamos:

1. Formular las hipótesis:

\(H_0: \beta_i=\beta_i^0\)

\(H_1:\beta_i \neq \beta_i^0\)

2. Calcular el estadístico de contraste:

\(t_i=\frac{\hat \beta_i - \beta_i}{S_{\hat \beta_i}}\)

3. Elegir el nivel de significación:

\(\alpha=0.05\) ó \(\alpha=0.1\)

4. Obtener en las tablas de la distribución t-student el valor crítico para un contraste bilateral (2 colas):

\(t_{n-k,0.025}\) ó \(t_{n-k,0.05}\)

Tambien se puede calcular p-valor asociado al valor de la t de student correspondiente al estadistico de contraste. El p-valor se define como la probabilidad de obtener un estadístico de contraste al menos tan extremo como el que evaluo.

5. Establecer la regla de decisión para aceptar la hipotesis nula \(H_0\):

Si \(\mid t_i \mid < t_{n-k,\frac{\alpha}{2}}\) acepto \(H_0\).

En el caso de que obtenga el p-valor del estadístico de contraste, este debe ser inferior al nivel de significación \(\alpha\) elegido.

Constraste de significación individual

Se denomina constraste de significación individual al contraste de hipótesis de que el parámetro \(\beta_i=0\). Es decir a:

1. Formular las hipótesis:

\(H_0: \beta_i=0\)

\(H_1:\beta_i \neq 0\)

2. Calcular el estadístico de contraste:

\(t_i=\frac{\hat \beta_i}{S_{\hat \beta_i}}\)

3. Elegir el nivel de significación:

\(\alpha=0.05\) ó \(\alpha=0.1\)

4. Obtener en las tablas de la distribución t-student el valor crítico o p-valor para un contraste bilateral (2 colas):

\(t_{n-k,0.025}\) ó \(t_{n-k,0.05}\)

5. Establecer la regla de decisión para aceptar la hipotesis nula \(H_0\):

Si \(\mid t_i \mid < t_{n-k,\frac{\alpha}{2}}\) acepto \(H_0\).

Aceptar \(H_0\) en este caso equivale a decir que el parámetro \(\beta_i\) es no significativo, ya que el valor cero entra dentro de los posibles valores que puede tomar el coeficiente en la población. En cambio rechazar \(H_0\), que ocurre cuando \(\mid t_i \mid > t_{n-k,\frac{\alpha}{2}}\) equivale a decir que el coeficiente \(\beta_i\) es significativo, y por tanto a adminitir algún tipo de influencia entre la variable \(X_i\) y la dependiente \(Y_i\).

Contraste de significación global(Tabla ANOVA)

La hipótesis de no significación global \(H_0:\beta_2=\beta_3=...=\beta_k=0\) se rechaza al nivel de significación \(\alpha\) construyendo el estadístico experimental: \[F_{exp}=\frac {\frac {SCE}{k-1}}{\frac {SCR}{n-k}}\]

y la regla de decisión que rechaza la hipótesis \(H_0\) ocurre cuando \(F_{exp}>F_{k-1,n-k,\alpha}\) .

El contraste de significación global se resume en el cuadro siguiente, en donde la variación total de la variable dependiente (\(SCT\)) se descompone en la explicada por la regresión (\(SCE\)) y en la no explicada (\(SCR\)). Los grados de libertad de estas tres sumas de cuadrados son \(n-1\),\(k-1\) y \(n-k\), respectivamente.

A partir de esta información muestral, podemos calcular el numerador y denominador del estadístico \(F_{exp}\).

Fuente de variación	Suma de cuadrados	Grados de libertad	Cuadrado medio	Estadistico \(F\)
Regresión	\(SCE\)	\(k-1\)	\(\frac {SCE}{k-1}\)	\(\frac {\frac {SCE}{k-1}}{\frac {SCR}{n-k}}\)
Residual	\(SCR\)	\(n-k\)	\(\frac {SCR}{n-k}\)
Total	\(SCT\)	\(n-1\)

Predicción

Una vez estimado y validado el modelo, una de sus aplicaciones más importantes consiste en poder realizar predicciones acerca del valor que tomaría la variable endógena en el futuro o para una unidad extramuestral. Esta predicción se puede realizar tanto para un valor individual como para un valor medio, o esperado, de la variable endógena, siendo posible efectuar una predicción puntual o por intervalos. Su cálculo se realiza mediante las expresiones que figuran a continuación:

1. Predicción individual: se trata de hallar el valor estimado para la variable Y un periodo hacia delante. En este caso basta con sustituir el valor de las variables exógenas en el modelo en el siguiente periodo y calcular el nuevo valor de Y:

\[\hat y_{t+1}= \hat \beta_1 x_{1,t+1} + ...+\hat \beta_k x_{k,t+1}=X_{t+1}\hat \beta\]

2. Intervalo de predicción. Para hallar un intervalo de predicción debe utilizarse la siguiente expresión:

\[\hat y_{t+1} \pm t_{n-k,\frac{\alpha}{2}}\hat \sigma \sqrt{1+X'_{t+1}(X'X)^{-1}X_{t+1}} \]

3. Intervalos de predicción para un valor medio o esperado,\(y_i\):

\[\hat y_{i} \pm t_{n-k,\frac{\alpha}{2}}\hat \sigma \sqrt{X'_{i}(X'X)^{-1}X_{i}} \]

Estimacion del modelo de regresión con R

R es un entorno especialmente diseñado para el tratamiento de datos, cálculo y desarrollo gráfico. Permite trabajar con facilidad con vectores y matrices y ofrece diversas herramientas para el análisis de datos.

R es una implementación open-source del lenguaje S (Bell Labs -principios de los 90), que también es la base del sistema S-Plus (entorno comercial). R y S-Plus aún comparten una gran mayoría de código e instrucciones, si bien R es software libre, gratuito en donde los usuarios disponen de libertad para ejecutar, copiar, distribuir, estudiar, cambiar y mejorar el software. De hecho R dispone de una comunidad de desarrolladores/usuarios detrás que se dedican constantemente a la mejora y a la ampliación de las funcionalidades y capacidades del programa. En la web http://www.r-project.org/ se encuentra disponible toda la información acerca de R. La instalación de R se realiza a través de la CRAN (ComprehensiveR Archive Network): http://cran.r-project.org

Actualmente R se distribuye para los siguientes Sistemas Operativos:

•Windows: entorno gráfico.

•Linux (Debian/Mandrake/SuSe/RedHat/VineLinux)

•MacOSX

•Código fuente: ampliación a sistemas Unix

Las funciones de R se agrupan en paquetes (packages, libraries), los que contienen las funciones más habituales se incluyen por defecto en la distribución de R, y el resto se encuentran disponibles en la Comprehensive R Archive Network (CRAN).

Las entidades que R crea y manipula se llaman objetos. Dichos objetos pueden ser:

•Escalares: números, caracteres, lógicos (booleanos), factores

•Vectores/matrices/listas de escalares

•Funciones

•Objetos ad-hoc

Dichos objetos se guardan en un workspace. Durante una sesión de R todos los objetos estarán en memoria, y se pueden guardar en disco para próximas sesiones.

R trabaja sobre estructuras de datos. La estructura más simple es un vector numérico, que consiste en un conjunto ordenado de números.

Un vector de reales se crea mediante la función c y se guarda con el nombre “Cantidad”.

Cantidad <-c(2.456,2.325,2.250,2.200,2.100,2.082,2.045,2.024)

Se crea ahora el vector de nombre “Precio”. mediana, se utilizan las siguientes funciones R:

Precio<-c(82,92,94,99,106,108,112,115)

Para obtener los estadísticos básicos del vector (Cantidad): media, desviación estandar, varianza y

mean(Cantidad)

## [1] 2.18525

sd(Cantidad)

## [1] 0.1515847

var(Cantidad)

## [1] 0.02297793

median(Cantidad)

## [1] 2.15

Si se quiere tener un resumen sumario de estadístico de una variable:

summary(Cantidad)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   2.024   2.073   2.150   2.185   2.269   2.456

En R los valores “desconocidos” o “no disponibles” (missings) se simbolizan con el valor especial NA (NotAvailable). Cualquier operación que incluya un NA en general devolverá NA como resultado.La función is.na nos permite saber si un elemento es missing o no.

Otros tipos de objectosen R.

•Arrays y matrices (matrix): generación multidimensional de los vectores. Todos los elementos de la matriz han de ser del mismo tipo.

•Factores (factor): útiles para el uso de datos categóricos.

•Listas (list): generalización de los vectores donde los elementos pueden ser de diferentes tipos (incluso vectores o nuevas listas).

•Data frames: matrices donde las diferentes columnas pueden tener valores de diferentes tipos.

•Funciones (function): conjunto de código de R ejecutable y parametrizable.

Una tabla debe estar en un objecto tipo matriz.

Ejemplo:

Tabla<-matrix(c(652,1537,598,242,36,46,38,21,218,327,106,67),nrow=3,byrow=T)

La función de R que nos permite estimar un modelo de regresión lineal es la función lm. La forma de invocar a la función para estimar un modelo de regresión lineal simple es lm(y~x). Se puede consultar la ayuda de la función para ver todas las posibilidades que ofrece.

En nuestro ejemplo, obtenemos:

lm(Cantidad~Precio)

## 
## Call:
## lm(formula = Cantidad ~ Precio)
## 
## Coefficients:
## (Intercept)       Precio  
##     3.53427     -0.01336

En lugar de invocar simplemente la función podemos guardar su resultado en una variable y veremos así que obtenemos más información.

reg = lm(Cantidad~Precio)

Si queremos obtener una gráfica con los resultados de la regresión realizada:

plot(Cantidad~Precio)
lines(reg$fitted~Precio)

Para realizar el análisis del modelo estimado utilizaremos la función summary. Así:

summary(reg)

## 
## Call:
## lm(formula = Cantidad ~ Precio)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.02875 -0.01359 -0.00154  0.01762  0.02574 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.5342726  0.0734707   48.10 5.41e-09 ***
## Precio      -0.0133567  0.0007235  -18.46 1.63e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02154 on 6 degrees of freedom
## Multiple R-squared:  0.9827, Adjusted R-squared:  0.9798 
## F-statistic: 340.8 on 1 and 6 DF,  p-value: 1.629e-06

Utilizando la base de datos “mtcars”, con datos sobre automoviles y sus características que incluye el programa R, cuya estructura se lista con la función “str”:

data(mtcars)
str(mtcars)

## 'data.frame':    32 obs. of  11 variables:
##  $ mpg : num  21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
##  $ cyl : num  6 6 4 6 8 6 8 4 4 6 ...
##  $ disp: num  160 160 108 258 360 ...
##  $ hp  : num  110 110 93 110 175 105 245 62 95 123 ...
##  $ drat: num  3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
##  $ wt  : num  2.62 2.88 2.32 3.21 3.44 ...
##  $ qsec: num  16.5 17 18.6 19.4 17 ...
##  $ vs  : num  0 0 1 1 0 1 0 1 1 1 ...
##  $ am  : num  1 1 1 0 0 0 0 0 0 0 ...
##  $ gear: num  4 4 4 3 3 3 3 4 4 4 ...
##  $ carb: num  4 4 1 1 2 1 4 2 2 4 ...

Realizamos una regresión lineal multiple entre el consumo de gasolina de cada coche y todas las caracteristicas que la base de datos incorpora para cada coche,

summary(lm(mpg~.,data=mtcars))

## 
## Call:
## lm(formula = mpg ~ ., data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4506 -1.6044 -0.1196  1.2193  4.6271 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 12.30337   18.71788   0.657   0.5181  
## cyl         -0.11144    1.04502  -0.107   0.9161  
## disp         0.01334    0.01786   0.747   0.4635  
## hp          -0.02148    0.02177  -0.987   0.3350  
## drat         0.78711    1.63537   0.481   0.6353  
## wt          -3.71530    1.89441  -1.961   0.0633 .
## qsec         0.82104    0.73084   1.123   0.2739  
## vs           0.31776    2.10451   0.151   0.8814  
## am           2.52023    2.05665   1.225   0.2340  
## gear         0.65541    1.49326   0.439   0.6652  
## carb        -0.19942    0.82875  -0.241   0.8122  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.65 on 21 degrees of freedom
## Multiple R-squared:  0.869,  Adjusted R-squared:  0.8066 
## F-statistic: 13.93 on 10 and 21 DF,  p-value: 3.793e-07

El sumario de estadística sobre la regresión nos ofrece los resultados de los contrastes de significación individual sobre los parámetros y el de significación global.

Para realizar una preducción del consumo de gasolina (mpg) en el modelo lineal que utiliza como explicativas el peso del vehículo (wt) y los caballos de vapor (hp), en un coche que pese, 3(1000lbs) y tenga 120 hp hay que utilizar la función R: “predict”

summary(lm(mpg~wt+hp,data=mtcars))

## 
## Call:
## lm(formula = mpg ~ wt + hp, data = mtcars)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.941 -1.600 -0.182  1.050  5.854 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 37.22727    1.59879  23.285  < 2e-16 ***
## wt          -3.87783    0.63273  -6.129 1.12e-06 ***
## hp          -0.03177    0.00903  -3.519  0.00145 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.593 on 29 degrees of freedom
## Multiple R-squared:  0.8268, Adjusted R-squared:  0.8148 
## F-statistic: 69.21 on 2 and 29 DF,  p-value: 9.109e-12

newdatos=data.frame(wt=3,drat=120)
predict(lm(mpg~wt+drat,data=mtcars),newdatos,interval="prediction",type="response")

##        fit       lwr      upr
## 1 189.0406 -158.0087 536.0898

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:

1. Estimar el modelo original y obtener la serie de residuos estimados

2. 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\]

3. 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.4493, df = 2, p-value = 0.2939

Autocorrelació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. Para facilitar la interpretación del test Durbin y Watson derivaron dos distribuciones: \(d_U\) y \(d_D\), 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 adopta la siguiente regla de decisión:

1. Si \(d \leq d_L\) rechazamos la hipótesis nula de no autocorrelación frente a la hipótesis alternativa de autocorrelación positiva.

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

3. Si \(d_U \leq d \leq (4-d_U)\) aceptamos la hipótesis nula de no autocorrelación.

4. 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 en donde 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.numeric

dwtest(y~x)

## 
##  Durbin-Watson test
## 
## data:  y ~ x
## DW = 1.6762, p-value = 0.09522
## alternative hypothesis: true autocorrelation is greater than 0

Deficiencias 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.50403957

Regresió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- numeric

plot(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.08372054

summary(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.06585964

MSPE.olsc

## [1] 0.06585964

MSPE.ridge

## [1] 0.04650589

Modelos con Variables Cualitativas

Modelos con Variables Cualitativas Explicativas

En un modelo econométrico, las variables representan a los conceptos u operaciones económicas que queremos analizar. Normalmente utilizamos variables cuantitativas, es decir, aquellas cuyos valores vienen expresados de forma numérica; sin embargo, también existe la posibilidad de incluir en el modelo econométrico información cualitativa, siempre que esta pueda expresarse de esa forma.

Las variables cualitativas expresan cualidades o atributos de los agentes o individuos (sexo, religión, nacionalidad, nivel de estudios, etc.) y también recogen acontecimientos extraordinarios como guerras, terremotos, climatologías adversas, huelgas, cambios políticos etc.

No cabe duda de que una forma de recoger factores de este tipo sería la utilización de variables proxy o aproximadas a las variables utilizadas. Por ejemplo, si quiero utilizar una variable que mida el nivel cultural de un país (variable cualitativa) puedo utilizar como variable proxy el número de bibliotecas existentes en un país, o representa una climatología adversa a partir de las temperaturas medias o precipitaciones. Sin embargo, no siempre es posible encontrar este tipo de variables y, en cualquier caso, debemos de ser conscientes de la posible existencia de errores en la definición de la variable.

Puesto que las variables cualitativas normalmente recogen aspectos de la presencia o no de determinado atributo (ser hombre o mujer, tener estudios universitarios o no tenerlos, etc.…) se utilizan variables construidas artificialmente, llamadas también ficticias o dummy, que generalmente toman dos valores, 1 ó 0, según se dé o no cierta cualidad o atributo. Habitualmente a la variable ficticia se le asigna el valor 1 en presencia de la cualidad y 0 en caso contrario. Las variables que toman valores 1 y 0, también reciben el nombre de variables dicotómicas o binarias.

Las variables dicotómicas pueden combinarse para caracterizar variables definidas por su pertenencia o no a un grupo. Si incluyo una variable cualitativa que me define la pertenencia o no de un país a un grupo, por ejemplo renta alta, media y baja, introduciré tres variables cualitativas en el modelo asociadas a la pertenencia o no a cada grupo; la primera caracterizaría a los individuos con renta alta, la segunda a los individuos con renta media, y la tercera a los individuos con renta baja.

Los modelos que utilizan variables cualitativas como regresores se diferencian en dos grupos, los modelos de Análisis de la Varianza o modelos ANOVA, que únicamente incluyen variables cualitativas como regresores; y los modelos de Análisis de la Covarianza o modelos ANCOVA que incluyen tanto variables cualitativas como cuantitativas. Los modelos ANOVA son muy utilizados en Sociología, Psicología, Educación, etc.; en Economía son más comunes los modelos ANCOVA.

Modelos ANOVA

Un problema estadístico clásico es la comparación de medias de dos distribuciones normales. Supongamos que las observaciones de la variable \(Y_i\) , provienen de dos distribuciones normales con medias \(\mu_1\) y \(\mu_2\) y varianza común \(\sigma^2\) . El tamaño de la primera distribución se circunscribe a las \(n_1\) primeras observaciones, y el de la segunda las \(n-n_1\) restantes observaciones. Queremos constrastar la hipótesis \(H_0:\mu_1=\mu_2\) frente a la alternativa \(H_0:\mu_1 \neq \mu_2\) al nivel de significación de \(\alpha\).

Este contraste de igualdad de medias cabe formularlo en el marco del modelo lineal general. Así, bajo \(H_0\) tenemos el siguiente modelo de regresión múltiple utilizando variables Dummy: \[ Y_i=\mu_1 D1_i + \mu_2 D2_i + u_i, (1)\]

donde \(D1_i\) toma valor uno en las \(n_1\) primeras observaciones y cero en las restantes, y \(D2_i\) toma cero en las \(n_1\) primeras observaciones y uno en las restantes, o lo que es lo mismo \(D2_i = 1-D1_i\).

Este modelo ANOVA puede formularse tambien a partir de las siguientes expresiones:

\[ Y_i=\beta_1 + \mu_1 D1_i + u_i, (2)\]

\[ Y_i=\beta_1 + \mu_2 D2_i + u_i,(3)\]

El contraste de igualdad de medias, se realiza a través del contraste de significación global, para el que construimos el estadístico experimental \(F_{exp}=\frac {\frac {R^2}{k-1}}{\frac {1-R^2}{n-k}}\), siendo el estadístico teórico \(F_{tco}\), la hipótesis se rechazaría con la regla de decisión \(F_{exp}>F_{tco}\).

Si se utiliza la especificación del modelo (2), el coeficiente asociado a la categoria \(D1\) dado por la suma (\(\beta_1+\mu_1\)), y para \(D2\) por \(\beta_1\). Si queremos contrastar la hipótesis de igualdad de medias en ambos grupos, equivaldría a contrastar la hipótesis nula \(H_0:\mu_1=0\).

Si se utiliza la especificación del modelo (3), el coeficiente asociado a la categoria \(D2\) dado por la suma (\(\beta_1+\mu_2\)), y para \(D2\) por \(\beta_1\). Si queremos analizar la hipótesis de ambos grupos, equivaldría a contrastar la hipótesis nula \(H_0:\mu_2=0\)

Partiendo de la base de datos “mtcars” preparamos un “Chunk” en el que construimos la tabla anova con la función “aov”, para los gastos por hogar y la variable categórica “am”, y estimamos un modelo ANOVA con la función “model.tables”.

##               Df Sum Sq Mean Sq F value   Pr(>F)    
## as.factor(am)  1  405.2   405.2   16.86 0.000285 ***
## Residuals     30  720.9    24.0                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Tables of means
## Grand mean
##          
## 20.09062 
## 
##  as.factor(am) 
##         0     1
##     17.15 24.39
## rep 19.00 13.00

La estimación MCO del modelo sería:

## 
## Call:
## lm(formula = mpg ~ 0 + as.factor(am), data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.3923 -3.0923 -0.2974  3.2439  9.5077 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## as.factor(am)0   17.147      1.125   15.25 1.13e-15 ***
## as.factor(am)1   24.392      1.360   17.94  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.902 on 30 degrees of freedom
## Multiple R-squared:  0.9487, Adjusted R-squared:  0.9452 
## F-statistic: 277.2 on 2 and 30 DF,  p-value: < 2.2e-16

## 
## Call:
## lm(formula = mpg ~ am, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.3923 -3.0923 -0.2974  3.2439  9.5077 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   17.147      1.125  15.247 1.13e-15 ***
## am             7.245      1.764   4.106 0.000285 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.902 on 30 degrees of freedom
## Multiple R-squared:  0.3598, Adjusted R-squared:  0.3385 
## F-statistic: 16.86 on 1 and 30 DF,  p-value: 0.000285

Modelos ANCOVA

Los modelos ANCOVA son modelos ANOVA que incluyen ademas de las variables cualitativas ó categóricas, variables cuantitativas entre las variables dependientes.

Partiendo de la especificación del modelo ANOVA del apartado anterior, al introducir la variable explicativa cuantitativa, \(X_i\), tendríamos:

\[ Y_i=\mu_1 D1_i + \mu_2 D2_i +\beta_2 X_i+ u_i\]

o sus equivalentes:

\[ Y_i=\beta_1 + \mu_1 D1_i +\beta_2 X_i+ u_i\]

\[ Y_i=\beta_1 + \mu_2 D2_i +\beta_2 X_i+ u_i\]

Para el análisis del comportamiento de cada grupo respecto a la pendiente, habría que especificar los siguientes modelos ANCOVA:

\[ Y_i=\delta_1 D1_iX_i + \delta_2 D2_iX_i + u_i\]

\[ Y_i=\beta_1 +\beta_2 X_i+ \delta_1 D1_i X_i+ u_i\]

\[ Y_i=\beta_1 +\beta_2 X_i+ \delta_2 D2_i X_i+ u_i\]

En los modelos ANOVA y ANCOVA se pueden incluir distintas variables categóricas. Consideremos una categorica \(E\) con tres categorías (1,2,3):

• \(E1_i\) que toma valor uno en los individuos pertenecientes a la categoría 1, y cero en los restantes

• \(E2_i\) que toma valor uno en los individuos pertenecientes a la categoría 2, y cero en los restantes

• \(E3_i\) que toma valor uno en los individuos pertenecientes a la categoría 3, y cero en los restantes.

Si bien a la hora de especificar el modelo hay que evitar los problemas de multicolinealidad entre todas las variables dummy incluidas y el término constante. Una forma de evitar los problemas es no incluir alguna de las categorías en forma de variable dummy, y dejar que la constante recoja el efecto de la categoría no incluida.

Una especificación posible de este modelo ANCOVA sería entonces:

\[ Y_i=\beta_1 +\beta_2 X_i+ \mu_2 D2_i + \alpha_1 E1_i + \alpha_2 E2_i + \alpha_3 E3_i u_i\]

Utilizando la base de datos “mtcars” vamos a estimar varios modelos ANCOVA:

## 
## Call:
## lm(formula = mpg ~ 0 + as.factor(am) + wt, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5295 -2.3619 -0.1317  1.4025  6.8782 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## as.factor(am)0  37.3216     3.0546  12.218 5.84e-13 ***
## as.factor(am)1  37.2979     2.0857  17.883  < 2e-16 ***
## wt              -5.3528     0.7882  -6.791 1.87e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.098 on 29 degrees of freedom
## Multiple R-squared:  0.9802, Adjusted R-squared:  0.9781 
## F-statistic: 478.1 on 3 and 29 DF,  p-value: < 2.2e-16

## 
## Call:
## lm(formula = mpg ~ am + wt, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5295 -2.3619 -0.1317  1.4025  6.8782 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 37.32155    3.05464  12.218 5.84e-13 ***
## am          -0.02362    1.54565  -0.015    0.988    
## wt          -5.35281    0.78824  -6.791 1.87e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.098 on 29 degrees of freedom
## Multiple R-squared:  0.7528, Adjusted R-squared:  0.7358 
## F-statistic: 44.17 on 2 and 29 DF,  p-value: 1.579e-09

## 
## Call:
## lm(formula = mpg ~ am + wt + wt * am, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.6004 -1.5446 -0.5325  0.9012  6.0909 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  31.4161     3.0201  10.402 4.00e-11 ***
## am           14.8784     4.2640   3.489  0.00162 ** 
## wt           -3.7859     0.7856  -4.819 4.55e-05 ***
## am:wt        -5.2984     1.4447  -3.667  0.00102 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.591 on 28 degrees of freedom
## Multiple R-squared:  0.833,  Adjusted R-squared:  0.8151 
## F-statistic: 46.57 on 3 and 28 DF,  p-value: 5.209e-11

## 
## Call:
## lm(formula = mpg ~ 0 + as.factor(am) + as.factor(gear) + wt, 
##     data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.5798 -2.4056 -0.3692  1.8198  5.7713 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## as.factor(am)0    35.0955     3.1862  11.015 1.72e-11 ***
## as.factor(am)1    35.2838     3.0091  11.726 4.20e-12 ***
## as.factor(gear)4   2.0769     1.7343   1.198    0.242    
## as.factor(gear)5  -1.0615     2.3845  -0.445    0.660    
## wt                -4.8782     0.7945  -6.140 1.46e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.968 on 27 degrees of freedom
## Multiple R-squared:  0.9831, Adjusted R-squared:  0.9799 
## F-statistic: 313.4 on 5 and 27 DF,  p-value: < 2.2e-16

La función R “lm”, presenta siempre una especificación que evita la multicolinealidad perfecta.

Realizar ejercicio libre con base de datos “npk”.

## 'data.frame':    24 obs. of  5 variables:
##  $ block: Factor w/ 6 levels "1","2","3","4",..: 1 1 1 1 2 2 2 2 3 3 ...
##  $ N    : Factor w/ 2 levels "0","1": 1 2 1 2 2 2 1 1 1 2 ...
##  $ P    : Factor w/ 2 levels "0","1": 2 2 1 1 1 2 1 2 2 2 ...
##  $ K    : Factor w/ 2 levels "0","1": 2 1 1 2 1 2 2 1 1 2 ...
##  $ yield: num  49.5 62.8 46.8 57 59.8 58.5 55.5 56 62.8 55.8 ...

Modelos con variables cualitativas endógenas.

Modelo probabilístico lineal

El modelo de probabilidad lineal se caracteriza por tener la variable endógena Y dicotómica o binaria, es decir toma el valor \(Y_i=1\) si un determinado suceso ocurre y el valor \(Y_i=0\) en caso contrario. Estos modelos están muy extendidos en el análisis estadístico pero encuentran una difícil aplicación en Economía debido a las dificultades de interpretación económica de los resultados que ofrecen este tipo de investigaciones. A este respecto, hay que considerar que estos modelos lo que realmente investigan es la probabilidad de que se dé una opción (valores \(Y_i=1\)) o no se dé (\(Y_i=0\)).

A pesar del carácter dicotómico de la variable endógena, el modelo de probabilidad lineal se especifica de la forma habitual, teniendo presente que las variables exógenas no son dicotómicas sino continuas:

\[ Y_i=\beta_1 + \beta_2 X_i + e_i, (4)\]

De acuerdo con la expresión (4), el hecho de que la variable endógena tome valores discretos (1 ó 0), el término de perturbación ei, puede tomar también dos valores únicamente:

• Si \(Y_i=0 \Rightarrow e_i=-\beta_1 - \beta_2X_i\) con probabilidad \(p\).

• Si \(Y_i=1 \Rightarrow e_i=1- \beta_1 - \beta_2X_i\) con probabilidad \((1-p)\).

Dado que la esperanza del término de error ha de ser nula \(E(e_i)=0\), entonces se demuestra que \(p=1-\beta_1 - \beta_2X_i\) y \((1-p) = \beta_1 + \beta_2X_i\), lo que permite evaluar la probabilidad de que la variable endógena tome el valor correspondiente:

• \(Prob (Y_i=0) = Prob (e_i=-\beta_1 - beta_2X_i=p=1-\beta_1 - \beta_2X_i\).

• \(Prob (Y_i=1) = Prob (e_i=1-\beta_1 - beta_2X_i=(1-p)=\beta_1 + \beta_2X_i\).

A su vez la varianza del término de perturbación, se calcularía a partir de:

\[Var(e_i)=p(1-p)=(1-\beta_1 - \beta_2X_i)(\beta_1 + \beta_2X_i)\]

Una problemática inherente a los estimadores MCO de estos modelos, es la siguiente:

• La perturbación aleatoria \(e_i\) no sigue una distribución Normal. Es sencillo observar este hecho ya que el carácter binario (1 ó 0) de la variable endógena afecta a la distribución de la perturbación, teniendo ésta una distribución Binomial . Este problema se atenúa cuando se utilizan tamaños de muestra (\(N\)) grandes en donde la distribución Binomial es susceptible de aproximarse a una Normal.

• La perturbación aleatoria no tiene una varianza constante (es heteroscedástica), lo cual supone una falta de eficiencia. Para solucionarlo habría que realizar transformaciones que nos diesen una perturbación homocedástica; esta transformación consiste en multiplicar todas las variables por una cierta cantidad que elimine el problema de la heteroscedasticidad. Dicha cantidad es:

\[\frac{1}{\sqrt{(1-\hat\beta_1 - \hat\beta_2X_i)(\hat\beta_1 + \hat\beta_2X_i)}}\]

siendo \(\hat\beta_1\) y \(\hat\beta_2\) los estimaciones MCO del modelo.

• No obstante, el mayor problema que plantean estos modelos es que las predicciones realizadas sobre la variable endógena no siempre se encuentran en el intervalo [0,1], ya que pueden ser mayores que cero y menores que uno. Este problema tiene dos soluciones, una es tomar como valor cero todas las estimaciones de la variable endógena con valores negativos, y uno cuando estas resulten mayores que uno; la segunda, solución es utilizar funciones de distribución que estén acotadas entre cero y uno como son la Logística y la Normal; de éstas se derivan los modelos Logit y Probit que pasamos a ver a continuación.

Modelo Logit

El problema que presentan los modelos probabilísticos lineales en cuanto a la existencia de predicciones establecidas fuera rango (negativas o mayores que uno), es debido a que utilizan una función de probabilidad que depende linealmente de las variables explicativas (\(X_i\)), que se resolverían acotando dicha distribución de probabilidad. El modelo Logit en concreto utiliza, para ello, la función de distribución logística:

Debido a que la función de distribución logística no tiene forma lineal, el modelo Logit se estima de forma diferente, así en vez de minimizar las sumas de las diferencias al cuadrado entre los valores observados y los estimados por el modelo, el carácter no lineal de los modelos Logit requiere la utilización del método de Máxima Verosimilitud para ser estimado, maximizando la verosimilitud de que un suceso tenga lugar, aunque se podría estimar por MCO mediante una transformación logarítmica de los datos (Gujarati, 1997).

La probabilidad de que \(Y_i=0\) (\(p\)) se define ahora mediante la siguiente expresión:

\[p=\frac{1}{1+e^{-z}}\]

donde \(z=\beta_1 + \beta_2X_i+e_i\).

La probabilidad de que \(Y_i=1\) (\(1-p\)) sería:

\[1-p=\frac{1}{1+e^{z}}\]

En consecuencia, la razón entre ambas será igual a:

\[\frac{p}{1-p}=\frac{1+e^{z}}{1+e^{-z}}=e^z\]

Tomando el logaritmo natural de la expresión anterior se obtiene:

\[L_i=\frac{p}{1-p}=loge^z=\beta_1 + \beta_2X_i\]

donde \(L_i\) es denominado Logit.

Los coeficientes \(\beta\) indican el cambio en el Logit causado por el cambio en una unidad en el valor de \(X_i\), mientras que los \(e^{\beta}\) definen el cambio en la razón de probabilidades \(\frac{p}{1-p}\) causado por el cambio en una unidad en el valor de \(X_i\). Si \(\beta\) es positivo,\(e^{\beta}\) será mayor que 1, es decir, \(\frac{p}{1-p}\) se incrementará;si \(\beta\) es negativo,\(e^{\beta}\) será menor que 1, y \(\frac{p}{1-p}\) disminuirá. Adicionalmente, puede demostrarse que el cambio en la probabilidad (\(p\)) causado por el cambio en una unidad en el valor de \(X_i\) es \(\beta(\frac{p}{1-p})\), es decir, depende no solo del coeficiente, sino también del nivel de probabilidad a partir del cual se mide el cambio.

A la hora de estimar un modelo Logit, hay que tener presente que para estimar el modelo además de los valores \(X_i\), se necesitan los valores del Logit (\(L_i\)), y por otro lado, la estimación de los coeficientes de modelo ha de realizarse utilizando el método de Máxima Verosimilitud, es decir, eligiendo como estimadores de los coeficientes  a aquellos que maximizan la función de verosimilitud, construida sobre la base de \(p=\frac{1}{1+e^{-z}}\). Pero si tenemos la posibilidad de agrupar los datos individuales, entonces podría estimarse el modelo por MCO.

Modelo Probit

Mientras que el modelo Logit utiliza la función de distribución logística para acotar la distribución de probabilidad en el modelo de probabilidad lineal, el modelo Probit utiliza la función de distribución Normal.

Las funciones de distribución normal y logística son muy semejantes: la diferencia principal es que la función de distribución normal se acerca más rápidamente a los ejes que la logística.

Para entender la filosofía del modelo Probit, vamos a suponer que existe una variable desconocida \(s\) que cumple lo siguiente:

Si \(I_i=\beta_1+\beta_2X_i \geq s\) entonces \(Y_i=1\)

Si \(I_i=\beta_1+\beta_2X_i<s\) entonces \(Y_i=0\)

Dado el supuesto de normalidad en un suceso, la probabilidad de que este sea menor o igual al valor (\(s\)), se calcula a partir de la función de distribución acumulada de una distribución Normal estandarizada, esto es, con esperanza cero y desviación típica uno.

\[pr(Y_i=1)=pr(I^*_i\geq I_i)=pr(\beta_1+\beta_2X_i \geq s)=F(\beta_1+\beta_2X_i)\]

\(F\) es la FDA normal estandar:

\[F(I_i)=\frac{1}{\sqrt{2\pi}}\int_{- \infty}^{\beta_1+\beta_2X_i}{e^{-\frac{z^2}{2}d(z)}}\]

donde \(Z_i\) es la variable normal estandarizada \(Z \sim N(0,\sigma^2)\)

Ahora para obtener información sobre \(I_i\), se toma la inversa de \(F(I_i)\):

\[I_i=F^{-1}(I_i)=F^{-1}(P)=\beta_1+\beta_2X_i\]

Al igual que ocurre en el probit si tenemos la posibilidad de agrupar los datos individuales, entonces podría estimarse el modelo por MCO.

Modelo Lineal Generalizado

Los modelos lineales (regresión, ANOVA, ANCOVA), se basan en los siguientes supuestos: 1. Los errores se distribuyen normalmente. 2. La varianza es constante. 3. La variable dependiente se relaciona linealmente con las variables independientes.

de manera analítica tendríamos: \[Y_i=\beta_1X_{1i}+\beta_2X_{1i}+...+\beta_kX_{1i}+u_i\], \(i=1,2,…, n\)

donde \(e_i\) esta distribuida de cómo una normal de media cero, varianza constante (homocedástica),y donde la covarianza entre \(e_i\) y \(e_j\) es nula para \(e_i\neq e_j\) (ausencia de autocorrelación). Es decir. Estos supuestos llevan implícito que la distribución de la variable dependiente \(Y_i\) sea también una normal \(Y_i\sim(\mu,\sigma^2)\).

donde \(\mu=\beta_1X_{1i}+\beta_2X_{1i}+...+\beta_kX_{1i}\)

En muchas ocasiones, sin embargo, nos encontramos con que uno o varios de estos supuestos no se cumplen por la naturaleza de la información. En algunos casos, estos problemas se pueden llegar a solucionar mediante la transformación de la variable respuesta (por ejemplo tomando logaritmos). Sin embargo estas transformaciones no siempre consiguen corregir la falta de normalidad, la heterocedasticidad (varianza no constante) o la no linealidad de nuestros datos.

Una alternativa a la transformación de la variable dependiente/respuesta y a la falta de normalidad es el uso de los modelos lineales generalizados (MLG).

Los MLG fueron formulados por John Nelder y Robert Wedderburn (1989) como una manera de unificar varios modelos estadísticos, incluyendo la regresión lineal, regresión logística y regresión de Poisson, bajo un solo marco teórico.

Los MLG son, por tanto, una extensión de los modelos lineales que permiten utilizar distribuciones no normales de los errores (binomiales, Poisson, gamma, etc) y varianzas no constantes.

Los MLG permiten especificar distintos tipos de distribución de errores, Cayuela (2010) expone los siguientes ejemplos:

• Poisson, muy útiles para conteos de acontecimientos, por ejemplo: número de heridos por accidentes de tráfico; número de hogares asegurados que dan parte de siniestro al día.

• Binomiales, de gran utilidad para proporciones y datos de presencia/ausencia, por ejemplo: tasas de mortalidad; tasas de infección; porcentaje de siniestros mortales.

• Gamma, muy útiles con datos que muestran un coeficiente de variación constante, esto es, en donde la varianza aumenta según aumenta la media de la muestra de manera constante, por ejemplo : número de heridos en función del número de siniestros

• Exponencial, muy útiles para los análisis de supervivencia.

Otra razón por la que un modelo lineal puede no ser adecuado para describir un fenómeno determinado es que la relación entre la variable respuesta y las variables independientes no es siempre lineal.

La función de vínculo se encarga de linealizar la relación entre la variable dependiente y las variables independientes mediante la transformación de la variable respuesta:

Tabla nº 1 Funciones de ligadura-vínculo más usadas

Función Vínculo	Formula	Uso
Identidad	\(\mu\)	Datos continuos con errores normales (regresión y ANOVA)
Logaritmica	\(log(\mu)\)	Conteos con errores de tipo de Poisson
Logit	\(log(\frac{\mu}{n-\mu})\)	Proporciones (datos entre O y 1) con errores binomiales
Recíproca	\(\frac{1}{\mu}\)	Datos continuos con errores Gamma
Cuadrada	\(\sqrt{\mu}\)	Conteos
Exponencial	\(\mu^n\)	Funciones de potencia

Fuente: (Cayuela, 2016)

Tabla nº2 Modelos MLG más comunes

Tipo de análisis	Variable respuesta	Variable explicativa	Función vínculo	Distribución de errores
Regresión	Continua	Continua	Identidad	Normal
ANOVA	Continua	Factor	Identidad	Normal
Regresión	Continua	Continua	Identidad	Gamma
Regresión	Conteo	Continua	Recíproca	Poisson
Tabla de contingencia	Conteo	Factor	Logarítmica	Poisson
Proporciones	Proporción	Continua	Logarítmica	Binomial
Regresión logística	Binaria	Continua	Logit	Binomial
Análisis de superviviencia	Tiempo	Recíproca	Identidad	Exponencial

Fuente: (Cayuela, 2016)

La estimación de los parámetros \(\beta\), se realiza por máximo verosimilitud,y los ajustes de \(\hat \mu_i\), se calculan como \(g^{-1}(x'_i\beta)\), una vez estimados los parámetros del vector \(\beta\). Para valorar el ajuste de los \(MLG\) se utiliza el estadístico chi-cuadrado, que se define como el doble de la diferencia entre el máximo del logaritmo de la verosimilitud que se podría conseguir con la mínima (o máxima) parametrización y el valor del máximo del logaritmo de la verosimilitud que se consigue con el modelo a evaluar, y el estadístico AIC (Akaike Information Criterion), formulado por Akaike (1974):

\[AIC=-2 \frac l N + 2 \frac k N \]

donde \(l\) es el valor en el óptimo del logaritmo de la función de verosimilitud con \(k\) parámetros estimados y \(N\) las observaciones. Siguiendo estos criterios, se seleccionará aquel modelo para el que se obtenga un AIC más bajo.

La especificación de un MLG se realiza en tres partes:

• La componente aleatoria correspondiente a la variable \(Y_i\) que sigue una distribución de la familia exponencial (normal, log-normal, poisson, gamma,…) • La componente sistemática, o predictor, que se denota \(\eta\) y corresponde al vector de \(n\) componentes \(\eta_i=\beta_1X_{1i}+\beta_2X_{1i}+...+\beta_kX_{1i}\) • La función de ligadura (o función link \(g(•)\)) que relaciona la esperanza matemática de la variable con el predictor lineal \(\eta_i=g(\mu_i)\).

Leemos base de datos de prestamos fallidos al consumo:

library(ISLR) 
str(Default) 

## 'data.frame':    10000 obs. of  4 variables:
##  $ default: Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 1 1 ...
##  $ student: Factor w/ 2 levels "No","Yes": 1 2 1 1 1 2 1 2 1 1 ...
##  $ balance: num  730 817 1074 529 786 ...
##  $ income : num  44362 12106 31767 35704 38463 ...

La variable dependiente será default, \(Y=YES\) e \(Y=No\).

Tenemos entonces tres variables predictivas: student, que caracteriza al consumidor como estudiante; balance, que es el saldo promedio de la tarjeta de crédito, e income que la renta del cliente.

Estimamos un modelo logit para explicar la admisión de alumnos, y evaluamos la tasa de

mylogit=glm(default~student+balance+income,family="binomial",data=Default) 
summary(mylogit) 

## 
## Call:
## glm(formula = default ~ student + balance + income, family = "binomial", 
##     data = Default)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4691  -0.1418  -0.0557  -0.0203   3.7383  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.087e+01  4.923e-01 -22.080  < 2e-16 ***
## studentYes  -6.468e-01  2.363e-01  -2.738  0.00619 ** 
## balance      5.737e-03  2.319e-04  24.738  < 2e-16 ***
## income       3.033e-06  8.203e-06   0.370  0.71152    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2920.6  on 9999  degrees of freedom
## Residual deviance: 1571.5  on 9996  degrees of freedom
## AIC: 1579.5
## 
## Number of Fisher Scoring iterations: 8

summary(mylogit$fitted.values)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0000103 0.0002798 0.0019662 0.0333000 0.0132244 0.9776263

fit.pred=ifelse(mylogit$fitted.values>0.5,1,0)
table(fit.pred,Default$default)

##         
## fit.pred   No  Yes
##        0 9627  228
##        1   40  105

Estimamos ahora un Probit para explicar la admisión de alumnos:

myprobit=glm(default~student+balance+income, data=Default,family=binomial(link=probit)) 
summary(myprobit) 

## 
## Call:
## glm(formula = default ~ student + balance + income, family = binomial(link = probit), 
##     data = Default)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.2226  -0.1354  -0.0321  -0.0044   4.1254  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -5.475e+00  2.385e-01 -22.960   <2e-16 ***
## studentYes  -2.960e-01  1.188e-01  -2.491   0.0127 *  
## balance      2.821e-03  1.139e-04  24.774   <2e-16 ***
## income       2.101e-06  4.121e-06   0.510   0.6101    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2920.6  on 9999  degrees of freedom
## Residual deviance: 1583.2  on 9996  degrees of freedom
## AIC: 1591.2
## 
## Number of Fisher Scoring iterations: 8

summary(myprobit$fitted.values)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0000000 0.0000183 0.0007723 0.0334850 0.0130713 0.9609657

fit.pred=ifelse(myprobit$fitted.values>0.5,1,0)
table(fit.pred,Default$default)

##         
## fit.pred   No  Yes
##        0 9639  238
##        1   28   95

Métodos de clasificación

Introducción

La clasificación supervisada es una de las tares que más frecuentemente son llevadas a cabo por los denominados Sistemas Inteligentes. Por lo tanto, un gran número de paradigmas desarrollados bien por la Estadística (Regresión Logística, Análisis Discriminante) o bien por la Inteligencia Artificial (Redes Neuronales, Inducción de Reglas, Árboles de Decisión, Redes Bayesianas) son capaces de realizar las tareas propias de la clasificación.

En el apartado anterior se han estudidado los métodos desarrollados por la estadística basados en el análisis de regresión: Regresión Logística y Probit, aquí estudiaremos otros métodos estadísticos como lo son, el Análisis Discriminantes y los K vecinos próximos, y los Arboles de Decisión, las Máquinas Soporte Vector, Redes Neuronales y el Clasificador Bayesiano desarrollados por la Inteligencia Artificial.

Paso previo a aplicar un método de clasificación, es la partición del conjunto de datos en dos conjuntos de datos más pequeños que serán utilizadas con los siguientes fines: entrenamiento y test . El subconjunto de datos de entrenamiento es utilizado para estimar los parámetros del modelo y el subconjunto de datos de test se emplea para comprobar el comportamiento del modelo estimado. Cada registro de la base de datos debe de aparecer en uno de los dos subconjuntos, y para dividir el conjunto de datos en ambos subconjuntos, se utiliza un procedimiento de muestreo: muestreo aleatorio simple o muestreo estratificado. Lo ideal es entrenar el modelo con un conjunto de datos independiente de los datos con los que realizamos el test.

Como resultado de aplicar un método de clasificación, se cometerán dos errores, en el caso de una variable binaria que toma valores 0 y 1, habrá ceros que se clasifiquen incorrectamente como unos y unos que se clasifiquen incorrectamente como ceros. A partir de este recuento se puede construir el siguiente cuadro de clasificación:

Valor real \(Y_i\)  Valor estimado \(\hat Y_i\)	\(Y_i=0\)	\(Y_i=1\)
\(\hat Y_i=0\)	\(P_{11}\)	\(P_{12}\)
\(\hat Y_i=1\)	\(P_{21}\)	\(P_{22}\)

Donde \(P_{11}\) y \(P_{22}\) corresponderán a predicciones correctas (valores 0 bien predichos en el primer caso y valores 1 bien predichos en el segundo caso), mientras que \(P_{12}\) y \(P_{21}\) corresponderán a predicciones erróneas (valores 1 mal predichos en el primer caso y valores 0 mal predichos en el segundo caso). A partir de estos valores se pueden definir los índices que aparecen en el siguiente cuadro:

Indice	Definición	Expresión
Tasa de aciertos	Cociente entre las predicciones correctas y el total de predicciones	\(\frac {P_{11}+P_{22}}{P_{11}+P_{12}+P_{21}+P_{22}}\)
Tasa de errores	Cociente entre las predicciones incorrectas y el total de predicciones	\(\frac {P_{12}+P_{21}}{P_{11}+P_{12}+P_{21}+P_{22}}\)
Especificidad	Proporción entre la frecuencia valores cero correctos y el total de valores cero observados	\(\frac {P_{11}}{P_{11}+P_{21}}\)
Sensibilidad	Proporción entre la frecuencia de valores uno correctos y el total de valores uno observados	\(\frac {P_{22}}{P_{12}+P_{22}}\)
Tasa de falsos ceros	Proporción entre la frecuencia de valores cero incorrectos y el total de valores cero observados	\(\frac {P_{21}}{P_{11}+P_{21}}\)
Tasa de falsos unos	Proporción entre la frecuencia de valores uno incorrectos y el total de valores uno observados	\(\frac {P_{12}}{P_{12}+P_{22}}\)

Un método para evaluar clasificadores alternativo a la métrica expuesta es la curva ROC (Receiver Operating Characteristic). La curva ROC es una representación gráfica del rendimiento del clasificador que muestra la distribución de las fracciones de verdaderos positivos y de falsos positivos. La fracción de verdaderos positivos se conoce como sensibilidad, sería la probabilidad de clasificar correctamente a un individuo cuyo estado real sea definido como positivo. La especificidad es la probabilidad de clasificar correctamente a un individuo cuyo estado real sea clasificado como negativo. Esto es igual a restar uno de la fracción de falsos positivos.

La curva ROC también es conocida como la representación de sensibilidad frente a (1-especificidad). Cada resultado de predicción representa un punto en el espacio ROC. El mejor método posible de predicción se situaría en un punto en la esquina superior izquierda, o coordenada (0,1) del espacio ROC, representando un 100% de sensibilidad (ningún falso negativo) y un 100% también de especificidad (ningún falso positivo). Una clasificación totalmente aleatoria daría un punto a lo largo de la línea diagonal, que se llama también línea de no-discriminación. En definitiva, se considera un modelo inútil, cuando la curva ROC recorre la diagonal positiva del gráfico. En tanto que en un test perfecto, la curva ROC recorre los bordes izquierdo y superior del gráfico. La curva ROC permite comparar modelos a través del área bajo su curva.

Curva ROC

En R existe una librería que ayuda a la representación de la curva ROC, el R-package ROCR.

library(ROCR)

## Loading required package: gplots

## 
## Attaching package: 'gplots'

## The following object is masked from 'package:stats':
## 
##     lowess

data(ROCR.simple)
pred <- prediction(ROCR.simple$predictions, ROCR.simple$labels)
perf <- performance(pred,"tpr","fpr")
plot(perf,colorize=TRUE)
abline(a=0, b= 1)

Validación cruzada

La validación cruzada o cross-validation es una técnica utilizada para evaluar los resultados de un análisis estadístico y garantizar que son independientes de la partición entre datos de entrenamiento y prueba. Consiste en repetir y calcular la media aritmética obtenida de las medidas de evaluación sobre diferentes particiones. Se utiliza en entornos donde el objetivo principal es la predicción y se quiere estimar cómo de preciso es un modelo que se llevará a cabo a la práctica.

Suponemos que tenemos un modelo con uno o más parámetros de ajuste desconocidos y unos datos de entrenamiento que queremos analizar. El proceso de ajuste optimiza los parámetros del modelo para que éste se ajuste a los datos de entrenamiento tan bien como pueda. Si tomamos una muestra independiente como dato de prueba (validación), del mismo grupo que los datos de entrenamiento, normalmente el modelo no se ajustará a los datos de prueba igual de bien que a los datos de entrenamiento. Esto se denomina sobreajuste y acostumbra a pasar cuando el tamaño de los datos de entrenamiento es pequeño o cuando el número de parámetros del modelo es grande.

Validación cruzada(De Joan.domenech91 - Trabajo propio, CC BY-SA 3.0, https://commons.wikimedia.org/w/index.php?curid=17617674)

Las maneras más habituales de realizar la validación cruzada son:

• Validación cruzada de K iteraciones o K-fold cross-validation. Los datos de muestra se dividen en K subconjuntos. Uno de los subconjuntos se utiliza como datos de prueba y el resto (K-1) como datos de entrenamiento. El proceso de validación cruzada es repetido durante k iteraciones, con cada uno de los posibles subconjuntos de datos de prueba. Finalmente se realiza la media aritmética de los resultados de cada iteración para obtener un único resultado. En la práctica, la elección del número de iteraciones depende de la medida del conjunto de datos. Lo más común es utilizar la validación cruzada de 10 iteraciones (10-fold cross-validation).

Validación cruzada dejando k interaciones (De Joan.domenech91 - Trabajo propio, CC BY-SA 3.0, https://commons.wikimedia.org/w/index.php?curid=17616792)

• Validación cruzada aleatoria.Este método consiste al dividir aleatoriamente el conjunto de datos de entrenamiento y el conjunto de datos de prueba. Para cada división la función de aproximación se ajusta a partir de los datos de entrenamiento y calcula los valores de salida para el conjunto de datos de prueba. El resultado final se corresponde a la media aritmética de los valores obtenidos para las diferentes divisiones. La ventaja de este método es que la división de datos entrenamiento-prueba no depende del número de iteraciones. Pero, en cambio, con este método hay algunas muestras que quedan sin evaluar y otras que se evalúan más de una vez, es decir, los subconjuntos de prueba y entrenamiento se pueden solapar.

Validación cruzada aleatoria (De Joan.domenech91 - Trabajo propio, CC BY-SA 3.0, https://commons.wikimedia.org/w/index.php?curid=17616794)

• Validación cruzada dejando uno fuera o Leave-one-out cross-validation (LOOCV). Se separan los datos de forma que para cada iteración tengamos una sola muestra para los datos de prueba y todo el resto conformando los datos de entrenamiento. La evaluación viene dada por el error, y en este tipo de validación cruzada el error es muy bajo, pero en cambio, a nivel computacional es muy costoso, puesto que se tienen que realizar un elevado número de iteraciones, tantas como N muestras tengamos y para cada una analizar los datos tanto de entrenamiento como de prueba.

Validación cruzada aleatoria (De Joan.domenech91 - Trabajo propio, CC BY-SA 3.0, https://commons.wikimedia.org/w/index.php?curid=17616792)

Análisis discriminante

El Análisis Discriminante (AD), introducido por Fisher (1936), es una técnica que se utiliza para predecir la pertenencia a un grupo (variable dependiente) a partir de un conjunto de predictores (variables independientes). El objetivo del AD es entender las diferencias de los grupos y predecir la verosimilitud de que una persona o un objeto pertenezca a una clase o grupo basándose en los valores que toma en los predictores. Ejemplos de análisis discriminante son distinguir entre innovadores y no innovadores de acuerdo a sus perfiles demográficos y sociales o el riesgo de impago de un préstamo a través de predictores económicos y sociodemográficos.

El análisis discriminante es conceptualmente muy similar al análisis de varianza multivariante de un factor. El AD trata de establecer una relación entre una variable dependiente no métrica (dicotómica o multidicotómica) y un conjunto de variables independientes métricas:

\[ Y_{1i}=X_{1i}+X_{2i}+...+X_{pi}\]

El propósito del AD consiste en aprovechar la información contenida en las variables independientes para crear una función Z combinación lineal de las p explicativas, capaz de diferenciar lo más posible a los k grupos. La combinación lineal para el análisis discriminante, función discriminante, se formula:

\[ Z_{jk}= \beta_0 + \beta_1 X_{1k} + \beta_2X_{2k}+...+ \beta_pX_{pk}\]

donde,

\(Z_{jk}\) es la puntuación \(Z\) discriminante \(j\) para el objeto \(k\).

\(\beta_0\) es el término constante.

\(\beta_i\) es la ponderación discriminante para la variable independiente \(i\).

\(X_{ik}\) es la variable independiente \(i\) para el objeto \(k\).

Una vez hallada la función discriminante, el resultado es una única puntuación Z discriminante compuesta para cada individuo en el análisis. Promediando las puntuaciones discriminantes para todos los individuos dentro de un grupo particular, obtenemos la media del grupo. Esta media es conocida como centroide. Cuando el análisis se realiza con dos grupos tenemos dos centroides, si es con tres serían tres los centroides, con k objetos tendremos k centroides.

En el caso de dos grupos y dos predictores o variables explicativas, la función discriminante es de la forma:

\[ Z_{jk}= \beta_0 + \beta_1 X_{1k} + \beta_2X_{2k}\]

Sustituyendo en la función discriminante el valor de las medias del grupo 1 en las variables \(X_1\) y \(X_2\), obtenemos el centroide del grupo 1:

\[ \bar Z_{1}= \beta_0 + \beta_1 \bar X_{11} + \beta_2 \bar X_{21}\]

De igual modo, sustituyendo las medias del grupo 2, obtenemos el centroide del grupo 2:

\[ \bar Z_{2}= \beta_0 + \beta_1 \bar X_{12} + \beta_2 \bar X_{22}\]

La función \(Z\) debe ser tal que la distancia entre los dos centroides sea máxima, consiguiendo de esta forma que los grupos estén lo más distantes posible. Podemos expresar esta distancia de la siguiente manera:

\[h= \bar Z_{1}-\bar Z_{2}\]

Es importante señalar que los grupos deben diferenciarse de antemano en las variables independientes. El análisis busca diferenciar los dos grupos al máximo combinando las variables independientes pero si los grupos no difieren en las variables independientes, no podrá encontrar una dimensión en la que los grupos difieran (Figura 5.1) .

Analisis Discriminante

Dicho de otro modo, si el solapamiento entre los casos de ambos grupos es excesivo, los centroides se encontrarán en la misma o parecida ubicación en el espacio p-dimensional y en esas condiciones, no será posible encontrar una función discriminante útil para la clasificación. Si los centroides están muy próximos, las medias de los grupos en la función discriminante serán tan parecidas que no será posible distinguir a los sujetos de uno y otro grupo.

La mayor utilidad de una función discriminante radica en su capacidad para clasificar nuevos casos. Ahora bien, la clasificación de casos es algo muy distinto de la estimación de la función discriminante. De hecho, una función perfectamente estimada puede tener una pobre capacidad clasificatoria.

Una vez obtenida la función discriminate podemos utilizarla, en primer lugar, para efectuar una clasificación de los mismos casos utilizados para obtener la función: esto permitirá comprobar el grado de eficacia la función desde el punto de vista de la clasificación. Si los resultados son satisfactorios, la función discriminante podrá utilizarse, en segundo lugar, para clasificar futuros casos de los que, conociendo su puntuación en las variables independientes, se desconozca el grupo al que pertenecen.

Una manera de clasificar los casos consiste en calcular la distancia existente entre los centroides de ambos grupos y situar un punto de corte \(z_0=\frac{\bar Z_{1}+\bar Z_{2}}{2}\) equidistante de ambos centroides. A partir de ese momento, los casos cuyas puntuaciones discriminantes sean mayores que el punto de corte,\(z_0\) serán asignados al grupo superior y los casos cuyas puntuaciones discriminantes sean menores que el punto de corte, \(z_0\), serán asignados al grupo inferior.

La regla de clasificación descrita sólo permite distinguir entre dos grupos, con lo que es difícilmente aplicable al caso de más de dos grupos e incluso a dos grupos con distinto tamaño, con tamaños desiguales es preferible utilizar una regla de clasificación que desplace el punto de corte hacia el centroide del grupo de menor tamaño buscando igualar los errores de clasificación. Para calcular este punto de corte se utiliza una distancia ponderada:

\(z_0=\frac{n_1\bar Z_{1}+n_2\bar Z_{2}}{n_1+n_2}\)

El AD solo admite variables cuantitativas como regresores, por lo que si alguna de las variables independientes es categórica, hay que utilizar otros métodos alternativos de clasificación.

La función R que realiza el Análisis Discriminante Lineal es “lda”. Para los 5 primeros datos, se dan los resultados de la clasificación (class), las probabilidades posteriores de pertenecer a la clase cero (posterior.0) o de pertenecer a la clase 1 (posterior.1), la probabilidad posterior es la probabilidad condicional que es asignada después de que la evidencia es tomada en cuenta.

Con la base de datos de R, Aid2, que contiene datos acerca de pacientes con sida,realizaremos un ejercicio de clasifición con AD, y evaluaremos los resultados con una métrica de porcentaje de aciertos y la curva ROC.

library(MASS)
data("Aids2")
str(Aids2)

## 'data.frame':    2843 obs. of  7 variables:
##  $ state  : Factor w/ 4 levels "NSW","Other",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ sex    : Factor w/ 2 levels "F","M": 2 2 2 2 2 2 2 2 2 2 ...
##  $ diag   : int  10905 11029 9551 9577 10015 9971 10746 10042 10464 10439 ...
##  $ death  : int  11081 11096 9983 9654 10290 10344 11135 11069 10956 10873 ...
##  $ status : Factor w/ 2 levels "A","D": 2 2 2 2 2 2 2 2 2 2 ...
##  $ T.categ: Factor w/ 8 levels "hs","hsid","id",..: 1 1 1 5 1 1 8 1 1 2 ...
##  $ age    : int  35 53 42 44 39 36 36 31 26 27 ...

# Lineal Discriminat Analisys
ad=lda(status~diag+age,data=Aids2)
#predicción
probs=predict(ad,newdata=Aids2,type="prob")
data.frame(probs)[1:5,]

##   class posterior.A posterior.D        LD1
## 1     A  0.54097063   0.4590294 -0.6423339
## 2     A  0.53868200   0.4613180 -0.6355380
## 3     D  0.03063147   0.9693685  2.0271192
## 4     D  0.03154984   0.9684502  2.0046295
## 5     D  0.09943774   0.9005623  1.1042254

table(probs$class,Aids2$status)

##    
##        A    D
##   A  794  332
##   D  288 1429

mean(probs$class==Aids2$status) #porcentaje de bien clasificados

## [1] 0.7819205

# Curva ROC
library(ROCR)
pred <- prediction(probs$posterior[,2],Aids2$status)
perf <- performance(pred,"tpr","fpr")
plot(perf,colorize=TRUE)
abline(a=0, b= 1)

Para realizar este ejercicio vamos a realizar una mineria de datos, dividiendo la base de datos en una muestra de entrenamiento y otra de test, la muestra de entrenamiento incluirá el 70% de las observaciones de la base de datos.

library(MASS)
data(iris)
attach(iris)
# división de la muestra en entrenamiento y validación
train=sample(seq(length(iris$Species)),length(iris$Species)*0.70,replace=FALSE)
# Lineal Discriminat Analisys
lda.tr=lda(Species[train]~.,data=iris[train,])
#predicción
probs=predict(lda.tr,newdata=iris[-train,],type="prob")
data.frame(probs)[1:5,]

##         class posterior.setosa posterior.versicolor posterior.virginica
## 6   virginica        0.3168932            0.3176922           0.3654145
## 13 versicolor        0.3473561            0.4328996           0.2197443
## 15 versicolor        0.2137919            0.4728890           0.3133191
## 18 versicolor        0.2796884            0.4174372           0.3028745
## 21 versicolor        0.3368992            0.4054189           0.2576819
##          x.LD1      x.LD2
## 6  -0.06303546  0.9674987
## 13 -0.96932719 -1.1596801
## 15 -1.46170745  1.5488047
## 18 -0.90995471  0.6303355
## 21 -0.75051528 -0.4973997

table(probs$class,iris$Species[-train])

##             
##              setosa versicolor virginica
##   setosa          1         12         8
##   versicolor      8          1         0
##   virginica       1          3        11

mean(probs$class==iris$Species[-train]) #porcentaje de bien clasificados

## [1] 0.2888889

Algorítmo K-vecinos más cercanos

El método K-nn (K nearest neighbors Fix y Hodges, 1951) es un método de clasificación supervisada (Aprendizaje, estimación basada en un conjunto de entrenamiento y prototipos) que sirve para estimar la función de densidad \(F(x / C_j)\) de las predictoraa \(x\) por cada clase \(C_j\).

Este es un método de clasificación no paramétrico, que estima el valor de la función de densidad de probabilidad o directamente la probabilidad a posteriori de que un elemento \(x\) pertenezca a la clase \(C_j\) a partir de la información proporcionada por el conjunto de prototipos o ejemplos. En el proceso de aprendizaje no se hace ninguna suposición acerca de la distribución de las variables predictoras.

K-vecinos

En la figura nº 5.2 se ilustra el funcionamiento de este método de clasificación. En la figura se encuentran representadas 12 muestras pertenecientes a dos clases distintas: la Clase 1 está formada por 6 cuadrados de color azul y la Clase 2 formada por 6 círculos de color rojo. En este ejemplo, se han seleccionado tres vecinos, es decir, (k=3).

De los 3 vecinos más cercanos a la muestra \(x\), representada en la figura por un aspa, uno de ellos pertenece a la Clase 1 y los otros dos a la Clase 2. Por tanto, la regla 3-nn asignará la muestra a la Clase 2. Es importante señalar que si se hubiese utilizado como regla de clasificación k=1, la 1-nn, la muestra sería asignada a la Clase 1, pues el vecino más cercano de la muestra pertenece a la Clase 1.

Un ejemplo de entrenamiento,\(x_i\), es un vector en un espacio característico multidimensional, que está descrito en términos de atributos, y pertenecerá a una de las clases de la clasificación. Los valores de los atributos del i-esimo ejemplo se representan por el vector -dimensional:

\[x_i=(x_{1i},x_{2i},...,x_{pi}) \epsilon X\]

El espacio es particionado en regiones por localizaciones y etiquetas de clases de los ejemplos de entrenamiento. Un punto en el espacio es asignado a la clase si esta es la clase más frecuente entre los k ejemplos de entrenamiento más cercano. Generalmente se usa la distancia euclidiana.

\[ d(x_i,x_j)=\sqrt{\sum_{r=1}^p(x_{ri}-x_{rj})^2} \]

La fase de entrenamiento del algoritmo consiste en almacenar los vectores característicos y las etiquetas de las clases de los ejemplos de entrenamiento. En la fase de test, la evaluación del ejemplo (del que no se conoce su clase) es representada por un vector en el espacio característico. Se calcula la distancia entre los vectores almacenados y el nuevo vector, y se seleccionan los k ejemplos más cercanos. El nuevo ejemplo es clasificado con la clase que más se repite en los vectores seleccionados.

El método k-nn supone que los vecinos más cercanos nos dan la mejor clasificación y esto se hace utilizando todos los atributos; el problema de dicha suposición es que es posible que se tengan muchos atributos irrelevantes que dominen sobre la clasificación, de manera que los atributos relevantes perderían peso entre otros veinte irrelevantes.

La mejor elección de k depende fundamentalmente de los datos; generalmente, valores grandes de k reducen el efecto de ruido en la clasificación, pero crean límites entre clases parecidas. Un buen k puede ser seleccionado mediante un procedimiento de optimización. El caso especial en que la clase es predicha para ser la clase más cercana al ejemplo de entrenamiento (cuando k=1) es llamada Nearest Neighbor Algorithm, Algoritmo del vecino más cercano.

Con la base de datos de R, Aid2, vamos a realizar una minería de datos, dividiendo la encuesta en dos muestras una de entrenamiento con el 70% de los datos y una muestra test con el 30% restante, se va a realizar el proceso con el primer vecino próximo k=1 (Nearest Neighbor Algorithm), para ello hay que instalar el package-R class, e invocar la función knn, dentro de esta librería la función knn permite elegir el numero de vecinos a aproximar, en está función todas las covariables han de ser numéricas, por este motivo seleccionamos en un “data.frame”" las variables numéricas de la base de datos. Evaluaremos los resultados con una métrica de porcentaje de aciertos.

library(class)
# Selección de variables
x=data.frame(Aids2[,3],Aids2[,7])
str(x)

## 'data.frame':    2843 obs. of  2 variables:
##  $ Aids2...3.: int  10905 11029 9551 9577 10015 9971 10746 10042 10464 10439 ...
##  $ Aids2...7.: int  35 53 42 44 39 36 36 31 26 27 ...

# Selección muestra entrenamiento
train=sample(seq(length(Aids2$status)),length(Aids2$status)*0.70,replace=FALSE)
# K-Nearest Neighbors
knn.prd=knn(x[train,],x[-train,],Aids2$status[train],k=3,prob=TRUE)
table(knn.prd,Aids2$status[-train])

##        
## knn.prd   A   D
##       A 198  91
##       D 128 436

mean(knn.prd==Aids2$status[-train]) #porcentaje de bien clasificados

## [1] 0.7432591

# Validación cruzada con la muestra de entrenamiento
knn.cross=knn.cv(x[train,],Aids2$status[train],k=3,prob=TRUE)
table(knn.cross,Aids2$status[train])

##          
## knn.cross    A    D
##         A  484  211
##         D  272 1023

Arboles de clasificación

Los árboles de decisión o clasificación tampoco son modelos estadísticos basados en la estimación de los parámetros de la ecuación propuesta, por tanto, no tenemos que estimar un modelo estadístico formal, son algoritmos para clasificar utilizando particiones sucesivas. Son apropiados cuando hay un número elevado de datos, siendo una de sus ventajas su carácter descriptivo que permite entender e interpretar fácilmente las decisiones tomadas por el modelo, revelando formas complejas en la estructura de datos que no se pueden detectar con los métodos convencionales de regresión.

Los árboles de decisión o de clasificación son un modelo surgido en el ámbito del aprendizaje automático (Machine Learning) y de la Inteligencia Artificial que partiendo de una base de datos, crea diagramas de construcciones lógicas que nos ayudan a resolver problemas. A esta técnica también se la denomina segmentación jerárquica. Es una técnica explicativa y descomposicional que utiliza un proceso de división secuencial, iterativo y descendente que partiendo de una variable dependiente, forma grupos homogéneos definidos específicamente mediante combinaciones de variables independientes en las que se incluyen la totalidad de los casos recogidos en la muestra.

Suponemos que se dispone de una muestra de entrenamiento que incluye la información del grupo al que pertenece cada caso y que sirve para construir el criterio de clasificación. Se comienza con un nodo inicial, dividiendo la variable dependiente a partir de una partición de una variable independendiente que se escoge de modo tal que de lugar a dos conjuntos homogeneos de datos (que maximizan la reducción en la impureza). Se elige, por ejemplo, la variable \(x_1\) y se determina un punto de corte, por ejemplo \(c\), de modo que se puedan separar los datos en dos conjuntos: aquellos con \(x_1 \leq c\) y los que tienen \(x_1 > c\). De este nodo inicial saldrán ahora dos: uno al que llegan las observaciones con \(x_1 \leq c\) y otro al que llegan las observaciones con \(x_1 > c\). En cada uno de estos nodos se vuelve a repetir el proceso de seleccionar una variable y un punto de corte para dividir la muestra. El proceso termina cuando se hayan clasificado todas las observaciones correctamente en su grupo.

En los árboles de decisión se encuentran los siguientes componentes: nodos, ramas y hojas. Los nodos son las variables de entrada, las ramas representan los posibles valores de las variables de entrada y las hojas son los posibles valores de la variable de salida. Como primer elemento de un árbol de decisión tenemos el nodo raíz que va a representar la variable de mayor relevancia en el proceso de clasificación. Todos los algoritmos de aprendizaje de los árboles de decisión obtienen modelos más o menos complejos y consistentes respecto a la evidencia, pero si los datos contienen incoherencias, el modelo se ajustará a estas incoherencias y perjudicará su comportamiento global en la predicción, es lo que se conoce como sobreajuste. Para solucionar este problema hay que limitar el crecimiento del árbol modificando los algoritmos de aprendizaje para conseguir modelos más generales. Es lo que se conoce como poda en los árboles de decisión.

Arboles de decisión

Las reglas de parada tratan de preguntar si merece la pena seguir o detener el proceso de crecimiento del árbol por la rama actual, se denominan reglas de prepoda ya que reducen el crecimiento y complejidad del árbol mientras se está construyendo:

• Pureza de nodo. Si el nodo solo contiene ejemplos o registros de una única clase se decide que la construcción del árbol ya ha finalizado.

• Cota de profundidad. Previamente a la construcción se fija una cota que nos marque la profundidad del árbol, cuando se alcanza se detiene el proceso.

• Umbral de soporte. Se especifica un número de ejemplos mínimo para los nodos, y cuando se encuentre un nodo con ejemplos por debajo del mínimo se para el proceso, ya que no consideramos fiable una clasificación abalada con menos de ese número mínimo de ejemplos.

Existen dos formas de poda muy comunes utilizadas en los diferentes algoritmos: la poda por coste-complejidad y la poda pesimista. En la poda por coste-complejidad se trata de equilibrar la precisión y el tamaño del árbol. La complejidad está determinada por el número de hojas que posee el árbol (nodos terminales). La poda pesimista utiliza los casos clasificados incorrectamente y obtiene un error de sustitución, eliminando los subárboles que no mejoran significativamente la precisión del clasificador.

Con la base de datos de R, iris, vamos a realizar una minería de datos, dividiendo la encuesta en dos muestras una de entrenamiento con el 70% de los datos y una muestra test con el 30% restante, se va ha realizar la clasificación utilizando arboles de decisión, para ello hay que instalar el package-R: “tree”, e invocar la función tree. Se realiza una poda por el procedimiento de coste-complejidad, y mediante un procedimiento de validación cruzada elegirá el mejor resultado. Para ello hay que invocar la función cv.tree con la opción FUN=prune.misclas. Evaluaremos los resultados con una métrica de porcentaje de aciertos. La salida ofrece la tasa de errores (Misclassification error rate).

library(tree)
# Selección muestra entrenamiento
train=sample(seq(length(iris$Species)),length(iris$Species)*0.70,replace=FALSE)
iris.tree = tree(iris$Species~.,iris,subset=train)
summary(iris.tree)

## 
## Classification tree:
## tree(formula = iris$Species ~ ., data = iris, subset = train)
## Variables actually used in tree construction:
## [1] "Petal.Length" "Sepal.Width" 
## Number of terminal nodes:  5 
## Residual mean deviance:  0.1041 = 10.41 / 100 
## Misclassification error rate: 0.01905 = 2 / 105

plot(iris.tree);text(iris.tree,pretty=0)

iris.tree

## node), split, n, deviance, yval, (yprob)
##       * denotes terminal node
## 
##  1) root 105 230.700 versicolor ( 0.32381 0.34286 0.33333 )  
##    2) Petal.Length < 2.45 34   0.000 setosa ( 1.00000 0.00000 0.00000 ) *
##    3) Petal.Length > 2.45 71  98.410 versicolor ( 0.00000 0.50704 0.49296 )  
##      6) Petal.Length < 4.85 36   9.139 versicolor ( 0.00000 0.97222 0.02778 )  
##       12) Petal.Length < 4.65 30   0.000 versicolor ( 0.00000 1.00000 0.00000 ) *
##       13) Petal.Length > 4.65 6   5.407 versicolor ( 0.00000 0.83333 0.16667 ) *
##      7) Petal.Length > 4.85 35   9.082 virginica ( 0.00000 0.02857 0.97143 )  
##       14) Sepal.Width < 2.65 5   5.004 virginica ( 0.00000 0.20000 0.80000 ) *
##       15) Sepal.Width > 2.65 30   0.000 virginica ( 0.00000 0.00000 1.00000 ) *

tree.pred=predict(iris.tree,iris[-train],type="class")
summary(tree.pred)

##     setosa versicolor  virginica 
##         50         49         51

with(iris[-train],table(tree.pred,Species))

##             Species
## tree.pred    setosa versicolor virginica
##   setosa         50          0         0
##   versicolor      0         46         3
##   virginica       0          4        47

# Mediante validación cruzada se busca el mejor arbol de decision
cv.iris=cv.tree(iris.tree,FUN=prune.misclass)
cv.iris

## $size
## [1] 5 3 2 1
## 
## $dev
## [1]  5  5 83 84
## 
## $k
## [1] -Inf    0   33   34
## 
## $method
## [1] "misclass"
## 
## attr(,"class")
## [1] "prune"         "tree.sequence"

plot(cv.iris)

El mismo ejercicio realizado con la función C5.0 del package C50.

library(C50)
data(iris)
 

# Selección de una submuestra de 105 (el 70% de los datos)
set.seed(101)
iris.indices <- sample(1:nrow(iris),size=105)
iris.entrenamiento <- iris[iris.indices,]
iris.test <- iris[-iris.indices,]
 
modelo <- C5.0(Species~., data = iris.entrenamiento,rules=TRUE)
modelo

## 
## Call:
## C5.0.formula(formula = Species ~ ., data = iris.entrenamiento, rules
##  = TRUE)
## 
## Rule-Based Model
## Number of samples: 105 
## Number of predictors: 4 
## 
## Number of Rules: 4 
## 
## Non-standard options: attempt to group attributes

# predicción
prediccion <- predict(modelo,newdata=iris.test)
# Matriz de confusión
tabla <- table(prediccion, iris.test$Species)
tabla

##             
## prediccion   setosa versicolor virginica
##   setosa         15          0         0
##   versicolor      0         11         0
##   virginica       0          1        18

# % correctamente clasificados
100 * sum(diag(tabla)) / sum(tabla)

## [1] 97.77778

La función, rpart, del package “rpart” tambien permite realizar arboles de decisión. Por defecto utiliza en metodo “Class”, que es el apropiado para estimar arboles de decisión. En la opción “control” encontramos parámetros opcionales para controlar el crecimiento del árbol. Por ejemplo, en control = rpart.control (minsplit = 30, cp = 0.001) se le indica que el número mínimo de observaciones en un nodo sea 30 antes de intentar una división y que una división debe pararse cuando el Coste factor de complejidad sea inferior a un 0.001.

En este caso utilizamos un ejemplo de niños que tienen problemas de espalda.

library(rpart)
fit <- rpart(Kyphosis ~ Age + Number + Start, data = kyphosis)
plot(fit)
text(fit, use.n = TRUE)

Realizar el ejercicio con la base de datos Aids2

Máquina Soporte Vector

Las Máquinas de Soporte Vectorial (Support Vector Machines SVMs) son un conjunto de algoritmos de aprendizaje supervisados que desarrollan métodos relacionados con los problemas de clasificación y regresión.

Como en la mayoría de los métodos de clasificación supervisada, los datos de entrada (los puntos) son vistos como un vector p-dimensional (una lista de p números). Dado un conjunto de puntos como un subconjunto de un conjunto mayor (espacio), en el que cada uno de ellos pertenece a una de dos posibles categorías, de manera que un algoritmo basado en SVM construye un modelo capaz de predecir si un punto nuevo (cuya categoría desconocemos) pertenece a una categoría o a la otra.

La SVM, intuitivamente, es un modelo que partiendo de un conjunto de ejemplos de entrenamiento, podemos etiquetarlos en diferentes clases y representar dichas muestras en puntos en el espacio para tratar de separar las diferentes clases mediante un espacio lo más amplio posible, para que cuando las nuevas muestras de los casos de test se pongan en correspondencia con dicho modelo puedan ser clasificadas correctamente en función de su proximidad.

En ese concepto de “separación óptima” es donde reside la característica fundamental de las SVM: este tipo de algoritmos buscan el hiperplano que tenga la máxima distancia (margen) con los puntos que estén más cerca de él mismo. Por eso también a veces se les conoce a las SVM como clasificadores de margen máximo. De esta forma, los puntos del vector que son etiquetados con una categoría estarán a un lado del hiperplano y los casos que se encuentren en la otra categoría estarán al otro lado.

La manera más simple de realizar la separación es mediante una línea recta, un plano recto o un hiperplano N-dimensional, pero los universos a clasificar no se suelen presentar en el ideal de las dos dimensiones como ocurre en el ejemplo gráfico, sino que un algoritmo SVM debe tratar con más de dos variables predictoras, curvas no lineales de separación, casos donde los conjuntos de datos no pueden ser completamente separados, clasificaciones en más de dos categorías. La representación por medio de funciones núcleo ó Kernel ofrece una solución a este problema, proyectando la información a un espacio de características de mayor dimensión el cual aumenta la capacidad computacional de la máquinas de aprendizaje lineal.

En la figura siguiente se presenta una programación R en donde SVM traza un hiperplano para clasificar la categoría de sertosa de la base de datos iris. Para ello hay que instalar el package-R: “e1017”, e invocar la función svm. Se estima un modelo con un lineal y un Kernel de base lineal (la función permite además funciones base polinomiales y sigmoides).

library("e1071")
library("rgl")
train <- iris
train$y <-ifelse(train[,5]=="setosa", 1, -1)
sv <- svm(y~Petal.Length+Petal.Width+Sepal.Length,data=train,kernel="linear",scale=FALSE,type="C-classification")
W <- rowSums(sapply(1:length(sv$coefs), function(i)sv$coefs[i]*sv$SV[i,]))
plot3d(train$Petal.Length, train$Petal.Width,train$Sepal.Length, col=ifelse(train$y==-1,"red","blue"),size = 2, type='s', alpha = .6)
rgl.planes(a = W[1], b=W[2], c=W[3],d=-sv$rho,color="gray", alpha=.4)

Clasificación de la clase Sertosa por SMV

Utilizamos la base de datos Aids, y realizamos un ejercicio de clasificación por svm, utilizando una muestra de entrenamiento y otra de test:

library(e1071)
# Selección de variables
datos=data.frame(Aids2[,-4])
str(datos)

## 'data.frame':    2843 obs. of  6 variables:
##  $ state  : Factor w/ 4 levels "NSW","Other",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ sex    : Factor w/ 2 levels "F","M": 2 2 2 2 2 2 2 2 2 2 ...
##  $ diag   : int  10905 11029 9551 9577 10015 9971 10746 10042 10464 10439 ...
##  $ status : Factor w/ 2 levels "A","D": 2 2 2 2 2 2 2 2 2 2 ...
##  $ T.categ: Factor w/ 8 levels "hs","hsid","id",..: 1 1 1 5 1 1 8 1 1 2 ...
##  $ age    : int  35 53 42 44 39 36 36 31 26 27 ...

# Selección muestra entrenamiento
train=sample(seq(length(Aids2$status)),length(Aids2$status)*0.70,replace=FALSE)
# se estima un modelo svm lineal para la muestra de entrenamiento
svmfit=svm(datos$status~.,data=datos,kernel="linear",scale=FALSE,subset=train)
print(svmfit)

## 
## Call:
## svm(formula = datos$status ~ ., data = datos, kernel = "linear", 
##     subset = train, scale = FALSE)
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  linear 
##        cost:  1 
##       gamma:  0.07142857 
## 
## Number of Support Vectors:  645

table(datos$status[train],svmfit$fitted)

##    
##       A   D
##   A 542 188
##   D 261 999

# Predicción para la muestra test
svm.pred=predict(svmfit,datos[-train,])
summary(svm.pred)

##   A   D 
## 375 478

with(datos[-train,],table(svm.pred,status))

##         status
## svm.pred   A   D
##        A 268 107
##        D  84 394

# se estima un modelo svm lineal para la muestra de entrenamiento y se predice la muestra de test
svmfit2=svm(datos$status~.,data=datos,kernel="radial",scale=FALSE,subset=train,probability=TRUE)
print(svmfit2)

## 
## Call:
## svm(formula = datos$status ~ ., data = datos, kernel = "radial", 
##     probability = TRUE, subset = train, scale = FALSE)
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  radial 
##        cost:  1 
##       gamma:  0.07142857 
## 
## Number of Support Vectors:  1913

svm.pred=predict(svmfit2,datos[-train,],probability=TRUE)
summary(svm.pred)

##   A   D 
## 191 662

with(datos[-train,],table(svm.pred,status))

##         status
## svm.pred   A   D
##        A 145  46
##        D 207 455

Red Neuronal Artificial

Una Red Neuronal Artificial (RNA) es un modelo matemático inspirado en el comportamiento biológico de las neuronas y en cómo se organizan formando la estructura del cerebro. Las redes neuronales intentan aprender mediante ensayos repetidos como organizarse mejor a si mismas para conseguir maximizar la predicción.

El primer modelo matemático de una neurona artificial, creado con el fin de llevar a cabo tareas simples, fue presentado en el año 1943 en un trabajo conjunto entre el psiquiatra y neuroanatomista Warren McCulloch y el matemático Walter Pitts. Un modelo de red neuronal se compone de nodos, que actúan como input, output o procesadores intermedios. Cada nodo se conecta con el siguiente conjunto de nodos mediante una serie de trayectorias ponderadas. Basado en un paradigma de aprendizaje, el modelo toma el primer caso, y toma inicial basada en las ponderaciones. Se evalúa el error de predicción y modifica las ponderaciones para mejorar la predicción, a continuación se evalúa el siguiente caso con las nuevas ponderaciones y se modifican para mejorar la predicción de los casos ya evaluados, el ciclo se repite para cada caso en lo que se denomina la fase de preparación o evaluación. Cuandos se ha calibrado el modelo, con la muesta test se evalúan los resultados globales.

Las redes neuronales se hicieron operativas por primera vez a finales de los 50. Rosenblat F. (1958) creó el perceptrón , un algoritmo de reconocimiento de patrones basado en una red de aprendizaje de computadora de dos capas usando una simple suma y la resta. A finales de los 60 el proceso de investigación sufrío un estancamiento, y ha sido a partir de los años 80 del siglo pasado, cuando se produjo el mayor desarrollo teórico. Son muchos los modelos de redes neuronales, el más utilizado es el algoritmo backpropagation fue creado por Werbos P. (1990).

Un ejemplo de modelo neuronal con \(n\) entradas, consta de:

• Un conjunto de entradas \(x_1,…x_n\).
• Los pesos sinápticos \(w_1,…w_n\), correspondientes a cada entrada.
• Una función de agregación,\(\sum\).
• Una función de activación,\(f\).
• Una salida,\(Y\).

Las entradas son el estímulo que la neurona artificial recibe del entorno que la rodea, y la salida es la respuesta a tal estímulo. La neurona puede adaptarse al medio circundante y aprender de él modificando el valor de sus pesos sinápticos, y por ello son conocidos como los parámetros libres del modelo, ya que pueden ser modificados y adaptados para realizar una tarea determinada.

En este modelo, la salida neuronal \(Y\) está dada por:

\(Y=f(\sum_{i=1}^{n}w_ix_i)\)

La función de activación se elige de acuerdo a la tarea realizada por la neurona. Entre las más comunes dentro del campo de las RNAs podemos destacar:

Funciones de Activación

Si qereemos un modelo de red neuronal para realizar tareas de clasificación en el plano (por lo que solo haremos uso de dos entradas, \(x_1\) y \(x_2\)). Para ello, vamos a considerar que existe una entrada de peso 1, llamada \(b\) y consideraremos como función de activación a la función signo definida por:

\(\Gamma(s) = \left\{ \begin{array}{cl} 1 & \mbox{si } s \geq 0 \\ -1 & \mbox{si } s < 0 \end{array} \right.\)

Por lo tanto, la salida neuronal \(Y\) estará dada en este caso por:

\(Y =\left\{ \begin{array}{cl} 1 & \mbox{si } w_1 x_1 + w_2 x_2 + b \geq 0 \\ -1 & \mbox{si } w_1 x_1 + w_2 x_2 + b < 0 \end{array} \right.\)

Supongamos que tenemos dos clases en el plano: la clase \(C_1\), formada por los círculos rojos, y la clase \(C_2\), formada por los círculos azules, donde cada elemento de estas clases está representado por un punto (x,y) en el plano. Supondremos además que tales clases son separables linealmente, es decir, es posible trazar una recta que separe estrictamente ambas clases.

Red Neuronal

Diremos que la neurona artificial clasifica correctamente las clases \(C_1\) y \(C_2\) si dados los pesos sinápticos \(w_1\) y \(w_2\) y el término aditivo \(b\), la recta con ecuación:

\(y = −\frac{w_1}{w_2} x − \frac{b}{w_2}\)

es una recta que separa las dos clases. La ecuación implícita de la recta es \(w_1x + w_2y + b = 0\).Obsérvese que si el punto \((x_0,y_0) \in C_1\) entonces \(w_1x_0+ w_2y_0+ b < 0\) y si \((x_0,y_0) \in C_2\) entonces \(w_1x_0+ w_2y_0+ b > 0\) Por lo tanto, dado el par \((x_0,y_0) \in C_1\cup C_2\), la neurona clasifica de la siguiente manera:

\((x_0,y_0) \in C_1 \Leftrightarrow Y = −1\)

\((x_0,y_0) \in C_2 \Leftrightarrow Y = 1\)

Si ahora tomamos dos clases \(C^∗_1\) y \(C^∗_2\) (separables linealmente) distintas a las anteriores, entonces la neurona puede no clasificar correctamente a estas clases, pues la recta anterior puede no ser una recta válida para separarlas. Sin embargo, es posible modificar los parámetros anteriores y obtener nuevos parámetros \(w^∗_1\),\(w^∗_2\) y \(b^∗\) tal que la recta \(y = −\frac{w^*_1}{w^*_2} x − \frac{b^*}{w^*_2}\), sirva ahora de separación entre ellos.

El proceso por el cual la neurona pasa de los parámetros \(w_1\),\(w_2\),\(b\) a los parámetros \(w^∗_1\),\(w^∗_2\), \(b^∗\) se conoce como método de aprendizaje. Este proceso es el que permite modificar los parámetros libres con el fin de que la neurona se adapte y sea capaz de realizar diversas tareas.

El método de aprendizaje que detallaremos a continuación y que utilizaremos para adaptar los parámetros libres con el fin de clasificar correctamente las clases \(C^∗_1\) y \(C^∗_2\) se conoce como método de error-corrección. Para aplicarlo es necesario:

• Un conjunto de entrenamiento, \(D = C_1\cup C_2\).
• Un instructor.
• Valores iniciales, \(w'_1\),\(w'_2\), \(b'\), arbitrarios de los parámetros libres.

El entrenamiento consiste en lo siguiente:

• El instructor toma un elemento \((x_0,y_0)\in D\) al azar y presenta éste a la neurona.
• Si la neurona clasifica mal este punto, es decir, si la salida de la neurona es \(Y=−1\) cuando \((x_0,y_0) \in C_2\) o \(Y=1\) cuando \((x_0,y_0) \in C_1\), entonces se aplica la siguiente corrección a los parámetros libres iniciales

\(w_1= w'_1+ d x_0\)

\(w_2= w'_2+ d y_0\)

\(b = b'+ d\)

donde el valor de \(d\) se obtiene de la siguiente manera:

\(d = =\left\{ \begin{array}{cl} 1 & \mbox{si } Y=-1 \mbox{ y } (x_0,y_0)\in C_2 \\ -1 & \mbox{si } Y=1 \mbox{ y } (X_0,y_0)\in C_1 \end{array} \right.\)

• Si la neurona clasifica bien el punto \((x_0,y_0)\), entonces no se realiza ninguna corrección.

• El procedimiento se repite pasando a la neurona otro punto del conjunto \(D\) y usando los últimos parámetros\(w_1\),\(w_2\), \(b\) corregidos (no los iniciales). Nuevamente, si la neurona clasifica mal el punto, entonces se aplica una corrección similar a la anterior.

• Esta tarea se repite con todos los puntos del conjunto \(D\).

• Si en el proceso hubo correcciones, entonces el procedimiento es repetido nuevamente con todos los puntos de \(D\).

• El entrenamiento termina cuando la neurona clasifica correctamente todos los elementos del conjunto de entrenamiento.

Un modelo neuronal utilizado para clasificación, cuya salida está dada de la forma anterior y que utiliza el método de error-corrección para modificar sus parámetros libres se conoce como Perceptron (el nombre deriva de la palabra en inglés “perception”). Estas neuronas pueden agruparse formando una RNA conocida como Perceptron múltiple.

Mas detalles en: http://www.cs.us.es/~fsancho/?e=72

Vamos a realizar una ejercicio de clasificación binaria con un conjunto de datos sobre fertilidad despues de abortos inducidos o expontaneos, que se encuentran el la librería “datasets”, la red neuronal que se utiliza por defecto es “rprot+”, relativa al algoritmo “Resilient Backpropagation”, que actualiza los pesos considerando únicamente el signo del cambio, es decir, si el cambio del error es en aumento (+) o disminución (-) entre una iteración y otra. Para detalles ver: https://en.wikipedia.org/wiki/Rprop

Los algoritmos de red que permite esta función R son:‘backprop’, ‘rprop+’, ‘rprop-’, ‘sag’, or ‘slr’.

La función para el calculo de error es la de entropía cruzada, que utiliza la distancia Kullback-Leibler.

library(neuralnet)

## 
## Attaching package: 'neuralnet'

## The following object is masked from 'package:ROCR':
## 
##     prediction

# Leemos los datos y hacemos dos conjuntos uno de entrenamiento y otro de test
data(infert, package="datasets")
str(infert)

## 'data.frame':    248 obs. of  8 variables:
##  $ education     : Factor w/ 3 levels "0-5yrs","6-11yrs",..: 1 1 1 1 2 2 2 2 2 2 ...
##  $ age           : num  26 42 39 34 35 36 23 32 21 28 ...
##  $ parity        : num  6 1 6 4 3 4 1 2 1 2 ...
##  $ induced       : num  1 1 2 2 1 2 0 0 0 0 ...
##  $ case          : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ spontaneous   : num  2 0 0 0 1 1 0 0 1 0 ...
##  $ stratum       : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ pooled.stratum: num  3 1 4 2 32 36 6 22 5 19 ...

n        <- nrow(infert)
muestra  <- sample(n,n*.70)
train    <- infert[muestra, ]
test     <- infert[-muestra, ]
# Entrenamos la red neuronal
net.infert <- neuralnet(case~parity+induced+spontaneous, train,err.fct="ce", linear.output=FALSE,likelihood=TRUE)
plot(net.infert)
# Computamos la red neuronal
resultados <- compute(net.infert,data.frame(test[,3:4],spontaneus=test[,6]))$net.result
# recodificamos los resultados y evaluamos
test.result=ifelse(resultados>0.5,1,0)
table(test.result,test$case)

##            
## test.result  0  1
##           0 38 14
##           1  5 18

Clasificador Bayesiano

Naïve Bayes es uno de los clasificadores más utilizados por su simplicidad y rapidez. Se trata de una técnica de clasificación y predicción supervisada que construye modelos que predicen la probabilidad de posibles resultados, en base al Teorema de Bayes, tambiéń conocido como teorema de la probalibidad condicionada:

\[P(A/B)=\frac{P(B/A)P(A)}{P(B)}\]

Un ejemplo nos servirá para entender la operativa de este clasificador. Supongamos un típico problema de clasificación en Minería de Datos, donde tenemos unos cuantos ejemplos que queremos clasificar. Cada ejemplo (también llamados instancias), se caracteriza por una serie de atributos. Además como estamos en clasificación supervisada, dichos ejemplos deben estar etiquetados, es decir, se debe conocer de qué clase son. Por tanto, tenemos un conjunto de 4 ejemplos con 3 atributos cada uno y 2 posibles clases.

Ejemplos	Atr. 1	Atr. 2	Atr. 3	Clase
\(x_1\)	1	2	1	positiva
\(x_2\)	2	2	2	positiva
\(x_3\)	1	1	2	negativa
\(x_4\)	2	1	2	negativa

Fuente:

Caben destacar dos partes en el algoritmo. La primera es la construcción del modelo, y la segunda clasificar un nuevo ejemplo con dicho modelo:

PARTE 1: Crear el modelo.

Para ello se necesitan cuatro pasos:

1. Calcular las probabilidades a priori de cada clase.

2. Para cada clase, realizar un recuento de los valores de atributos que toma cada ejemplo. Se debe distribuir cada clase por separado para mayor comodidad y eficiencia del algoritmo.

3. Aplicar la Corrección de Laplace, para que los valores “cero” no den problemas.

4. Normalizar para obtener un rango de valores [0,1].

A continuación se detalla cada paso:

Paso 1: Hay varias formas de establecer la probabilidad a priori a cada clase. La más intuitiva es que todas ellas tengan la misma probabilidad, es decir, 1 divido entre el número de clases. Si contamos con la opinión de un experto, puede que éste nos proporciones dichas probabilidades.

En este caso, vamos a considerar cada clase equiprobable. Por tanto:

\(P(positivo) = 1/2 = 0,5\)

\(P(negativo) = 1/2 = 0,5\)

Paso 2: Se crea una tabla (un otro contenedor) para cada clase y se hace un recuento de valores:

positivo	valor 1	valor 2
atr.1	1	1
atr.2	0	2
atr.3	1	1

negativo	valor 1	valor 2
atr.1	1	1
atr.2	2	0
atr.3	0	2

Paso 3: Se suma una unidad a cada valor de la tabla anterior.

positivo	valor 1	valor 2
atr.1	2	2
atr.2	1	3
atr.3	2	2

negativo	valor 1	valor 2
atr.1	2	2
atr.2	3	1
atr.3	1	3

Paso 4: Se normalizan todos los valores de la tabla del siguiente modo. Cada celda se divide por la suma de los valores de la fila. Por ejemplo, el valor: “valor 1” del atributo: “atr. 2” se actualiza con: 1 / (1+3) = 0,25.

positivo	valor 1	valor 2
atr.1	0.5	0.5
atr.2	0.25	0.75
atr.3	0.5	0.5

negativo	valor 1	valor 2
atr.1	0.5	0.5
atr.2	0.75	0.25
atr.3	0.25	0.75

PARTE 2: Clasificar un nuevo ejemplo. Por ejemplo, \(x_5 = (1,1,1)\).

Para ello se necesitan dos pasos:

1. Para cada clase disponible, se determinan los valores de probabilidad de cada valor de los atributos del nuevo ejemplo.

2. Aplicar la fórmula de Naïve Bayes.

Si generalizamos el teorema de Bayes a \(n\) atributos, obtenemos la siguiente expresión:

\(P(A | b_1,b_2,...,b_n)=P(A) P(b_1| A) P(b_2| A),...,P(b_n | A) = P(A)P(b_j| A)\)

que resolvemos:

\(Solucion=argmax_{i=1}^n P(A) P(b_i| A)\)

A continuación se detalla cada paso:

Paso 1: En este caso, se para cada tabla escogemos los valores: 1 de cada atributo, es decir, la primera columna de cada tabla, ya que todos los valores del ejemplo \(x_5\) son: \(1\).

positivo	valor 1	valor 2
atr.1	0.5	0.5
atr.2	0.25	0.75
atr.3	0.5	0.5

negativo	valor 1	valor 2
atr.1	0.5	0.5
atr.2	0.75	0.25
atr.3	0.25	0.75

Paso 2: Se aplica la fórmula:

\(P(positivo | x_5 ) = 0,5 • ( 0,5 • 0,25 • 0,5 ) = 0,03125\) \(P( negativo | x_5 ) = 0,5 • ( 0,5 • 0,75 • 0,25 ) = 0,046875\) Como la clase más probable es “negativa”, clasificamos el nuevo ejemplo (\(x_5\)) como negativo.

Esta forma de clasificar, se denomina Tecnica de maximo a posteriori.

A continuación se presenta un ejemplo de clasificar clientes con el método Naive Bayes, para clasificar a los clientes según su probabilidad de compra (Si o No) y en base a determinados atributos personales(Estado Cvivil,Profesion, etc.).

El cálculo de la probabilidad de cada atributo y para cada cliente nuevo seria algo como:

La programación en R sería la siguiente:

library(e1071)
data(iris)
m <- naiveBayes(Species ~ ., data = iris)
## alternatively:
m <- naiveBayes(iris[,-5], iris[,5])
m

## 
## Naive Bayes Classifier for Discrete Predictors
## 
## Call:
## naiveBayes.default(x = iris[, -5], y = iris[, 5])
## 
## A-priori probabilities:
## iris[, 5]
##       setosa   versicolor    virginica 
## 0.3333333333 0.3333333333 0.3333333333 
## 
## Conditional probabilities:
##             Sepal.Length
## iris[, 5]     [,1]         [,2]
##   setosa     5.006 0.3524896872
##   versicolor 5.936 0.5161711471
##   virginica  6.588 0.6358795933
## 
##             Sepal.Width
## iris[, 5]     [,1]         [,2]
##   setosa     3.428 0.3790643691
##   versicolor 2.770 0.3137983234
##   virginica  2.974 0.3224966382
## 
##             Petal.Length
## iris[, 5]     [,1]         [,2]
##   setosa     1.462 0.1736639965
##   versicolor 4.260 0.4699109772
##   virginica  5.552 0.5518946957
## 
##             Petal.Width
## iris[, 5]     [,1]         [,2]
##   setosa     0.246 0.1053855894
##   versicolor 1.326 0.1977526800
##   virginica  2.026 0.2746500556

table(predict(m, iris), iris[,5])

##             
##              setosa versicolor virginica
##   setosa         50          0         0
##   versicolor      0         47         3
##   virginica       0          3        47

Combinado clasificadores

Recientemente en el área de la Inteligencia Artificial el concepto de combinación de clasificadores ha sido propuesto como una nueva dirección para mejorar el rendimiento de los clasificadores individuales. Estos clasificadores pueden estar basados en una variedad de metodologías de clasificación, y pueden alcanzar diferentes ratios de individuos bien clasificados. El objetivo de la combinación de clasificadores individuales es el ser más certeros, precisos y exactos. Los métodos multiclasificadores más conocidos son el Bagging (Breiman, 1966) y Boosting (Freund y Schapire, 1996).

El método propuesto por Breinan (1966) intenta aunar las características del Boostrapping y la agregación incorporando los beneficios de ambos (Boostrap AGGregatiNG). La operativa del método es la siguiente:

• Se generan muestras aleatorias que serán los conjuntos de entrenamiento. Las muestras se generan a través de un muestreo aleatorio con reemplazamiento.

• Cada subconjunto de entrenamiento aprende un modelo.

• Para clasificar un ejemplo se predice la clase de ese ejemplo para cada clasificador y se clasifica en la clase con mayor voto.

El método propuesto por Freund y Schapire (1996), está basado en la asignación de un peso a cada conjunto de entrenamiento. Cada vez que se itera se aprende un modelo que minimiza la suma de los pesos de aquellos ejemplos clasificados erróneamente. Los errores de cada iteración sirven para actualizar los pesos del conjunto de entrenamiento, incrementando el peso de los mal clasificados y reduciendo el peso de aquellos que han sido correctamente clasficados. La decisión final para un nuevo patrón de clasificación viene dada por la votación mayoritaria ponderada entre los diferentes conjuntos de entrenamiento.

El package R: “ipred” opera multiclasificadores por los métodos bagging y boosting.

Bagging (o Bootstrap Aggregating), consiste en crear diferentes modelos usando muestras aleatorias con reemplazo y luego combinar o ensamblar los resultados.

La técnica de Bagging sigue estos pasos:

1. Divide el set de Entrenamiento en distintos sub set de datos, obteniendo como resultado diferentes muestras aleatorias con las siguientes características:

• Muestra uniforme (misma cantidad de individuos en cada set)

• Muestras con reemplazo (los individuos pueden repetirse en el mismo set de datos).

• El tamaño de la muestra es igual al tamaño del set de entrenamiento, pero no contiene a todos los individuos ya que algunos se repiten.

• Si se usan muestras sin reemplazo, suele elegirse el 50% de los datos como tamaño de muestra

2. Luego se crea un modelo predictivo con cada set, obteniendo modelos diferentes

3. Luego se construye o ensambla un único modelo predictivo, que es el promedio de todos los modelos.

Red Neuronal Artificial Precio Vivienda

Sobre el Bagging hay que tener presente (http://apuntes-r.blogspot.com.es/search/label/Bagging):

• Disminuye la varianza de un data set al realizar remuestreo con reemplazo..

• Si no existe varianza en el data set, la tecnica de Baggin no mejora significativamente el modelo.

• Es recomendable en modelos de alta inestabilidad (data set con mucha varianza). Ejemplo de inestabilidad: el % de error de la predicción de fraudes de enero , es muy diferente al de febrero.

• Mientras más inestable es un modelo, mejor será la predicción al usar Bagging.

• Se reduce el overfetting o sobre entrenamiento de modelos. Esto porque los modelos no pueden sobreaprender o memorizar ya que ninguno tiene todos los datos de entrenamiento.

• Mejora la predicción, ya que lo que no detecta un modelo lo detectan los otros.

• Reduce el ruido de los outliers, ya los outliers no pueden estar presenten en todos los modelos.

• No mejora significativamente las funciones lineales, ya que el ensamble de una función lineal da como resultado otra función linear.

• Una técnica mejorada del Bagging es el Random Forest que ademas de elegir un grupo aleatorio de individuos, también elige un grupo aleatorio de variables.

• Los diferentes modelos creados con la técnica Bagging pueden considerarse como algoritmos que buscan respuestas (o hipótesis) en un data set (o espacio h). Como cada algoritmo tiene un set de datos diferentes, cada uno creará una hipótesis diferentes sobre la realidad.

Realizamos un bagging para un clasificador binario, utilizando un scrip obtendido en: http://apuntes-r.blogspot.com.es/2015/01/bagging-para-clasificador-binario.html#more

# Bagging para un clasificador binario
 
 
# -----------------------------------------------------------
# PASO 1. Carga de package y datos
suppressMessages(library(C50));data(churn);   # para datos
suppressMessages(library(foreach))            # para la logica de bagging
library(rpart)                                # para algoritmo Arbol Decision
 
Train        <- churnTrain[,1:20]
Test         <- churnTest [,1:20]
Clase        <- unique(churnTrain$churn)
Iteraciones  <- 5  # Debe ser impar, para evitar el empate en el "Voto mayoritario"
 
#-----------------------------------------------------------
# PASO 2. Crea modelos y crea predicciones en columnas
Prediccion   <- foreach(i=1:Iteraciones,.combine=cbind) %do% {
muestra    <- sample(nrow(Train), size=floor((nrow(Train)*.7)))
modelo     <- rpart(churn ~., data=Train[muestra,])
as.character(predict(modelo, Test, type="class"))
}  
 
# -----------------------------------------------------------
# PASO 3. Evaluación del Voto Mayoritario
Prediccion                 <- as.data.frame(Prediccion)
Prediccion$Cantidad_yes    <- rowSums(Prediccion [, 1:Iteraciones] == "yes")
Prediccion$Cantidad_no     <- rowSums(Prediccion [, 1:Iteraciones] == "no")
Prediccion$Prediccion      <- with(Prediccion,ifelse(Cantidad_yes>Cantidad_no,"yes","no"))
MC                         <- table(Test[, "churn"],Prediccion[,"Prediccion"])
Efectividad                <- MC["yes","yes"]/(MC["yes","yes"]+MC["yes","no"])
print(MC)

##      
##         no  yes
##   yes   82  142
##   no  1430   13

print(Efectividad)

## [1] 0.6339285714

Agrupación de la información

Análisis de componentes principales

El Análisis de Componentes Principales (ACP) es una técnica estadística de síntesis de la información, o reducción de la dimensión (número de variables.

Dadas \(n\) observaciones de \(p\) variables, se buscan \(m<p\) variables que sean combinaciones lineales de las \(p\) originales y que estén incorreladas, recogiendo la mayor parte de la información o variabilidad de los datos. Hay que tener presente que, si las variables originales están incorreladas de partida, entonces no tiene sentido realizar un análisis de componentes principales.

Se considera una serie de variables \((x_1,x_2,...x_p)\) sobre un grupo de objetos o individuos y se trata de calcular, a partir de ellas, un nuevo conjunto de variables \((y_1,y_2,...y_p)\),incorreladas entre sí, cuyas varianzas vayan decreciendo progresivamente. Cada \(y_j\) (donde \(j=1,...,p\)) es una combinación lineal de las \(x_1,x_2,...x_p\) originales, es decir:

\(y_j=a_{j1}x_1+a_{j2}x_2+...+a_{jP}x_p=a'_jx\)

donde \(a'_j=(a_{j1},a_{j2},...,a_{jP})\) y

\[x= \begin{bmatrix} x_{1}\\ x_{2}\\ . . . x_{p} \end{bmatrix}\]

El objeto del análisis es maximizar la varianza del componente obtenido, pero manteniendo la ortogonalidad de la transformación, esto es \(a'_ja_j=\sum_{j=1}^p a_{kj}^2=1\).

El primer componente por tanto se obtendrá:

\(Var(y_1)=Var(a'_1 x)=a'_1\sigma a_1\)

El problema consiste en maximizar la función \(a'_1\sigma a_1\) donde \(\sigma\) es la matriz de covarianzas de las \(p\) variables,sujeta a la restricción \(a'_1a_1=1\), utilizandose para ello los multiplicadores de Lagrange:

\(L(a_1)=a'_1\sigma a_1-\lambda (a'_1a_1-1)\)

siendo ahora \(a_1\) la incognita, se busca el maximo:

\(\frac{L(\partial a_1)}{\partial a_1}=0 \longrightarrow (\sigma - \lambda I)a_1=0\)

La solución del sistema de ecuaciones requiere que el determinante \(\mid \sigma - \lambda I \mid\) sea no nulo, de manera que \(\lambda\) es un autovalor de \(\sigma\), que al ser definida positiva, tendrá \(p\) autovalores distintos.

Entonces, para maximizar la varianza de \(y_1\) se tiene que tomar el mayor autovalor, \(\lambda_1\) y el correspondiente autovector de \(a_1\).

El segundo componente principal,\(y_2=a'_2x\) ha de estar incorrelacionado con \(y_1\), es decir:

\(cov(y_1,y_2)=cov(a'_1x,a'_2x)=a'_1 \sigma a_2=0\).

Equivale a:

\(a'_1 \lambda a_2=0\)

Y a maximizar las \(Var(y_2)=a'_2\sigma a_2\) imponinedo ahora dos restricciones:

\(a'_2a_2=1\)

\(a'_2a_1=0\)

Se forma el lagrangiano \(L(a_2)=a'_12\sigma a_2-\lambda (a'_2a_2-1) - \delta (a'2a_1)\), se deriva respecto a \(a_2\) y se iguala a cero:

\(\frac{L(\partial a_2)}{\partial a_2}=0 \longrightarrow (\sigma - \lambda I)a_2=0\).

Elegimos \(\lambda_2\) como segundo autovalor de la matriz \(\sigma\) asociado al autovector asociado \(a_2\).

Entonces todos los p componentes de \(y\) se pueden expresar como el producto de una matriz formada por los autovectores, multiplicada por el vector \(x\) que contiene las variables originales:

\(y=Ax\)

donde la \(Var(y_j)=\lambda_j\), y la \(Cov(y_j,y_s)=0\).

Si sumamos todos los autovalores,tendremos entonces la varianza total de los componentes, que por las propiedades del operador traza, acaba siendo igual que la varianza de los variables originales:

\(\sum Var(y_j)= \sum (\lambda_j)=\sum Var(x_j)\)

El porcentaje de varianza total que recoge cada componentes principales, es el critierio utilizado para reducir las \(p\) variables de nuestro conjunto de datos, por \(m\) combinaciones lineales de dichas variables.

El ACP se emplea sobre todo en análisis exploratorio de datos (Analisis factorial) y para construir modelos predictivos cuando las variables explicativas tienen problemas de colinealidad.

El ACP comporta el cálculo de la descomposición en autovalores de la matriz de covarianza, normalmente tras centrar los datos en la media de cada atributo.

El análisis de componentes principales (PCA) puede realizarse con una función de la librería básica de «R»: prcomp o princomp.

require(graphics)
str(USArrests)

## 'data.frame':    50 obs. of  4 variables:
##  $ Murder  : num  13.2 10 8.1 8.8 9 7.9 3.3 5.9 15.4 17.4 ...
##  $ Assault : int  236 263 294 190 276 204 110 238 335 211 ...
##  $ UrbanPop: int  58 48 80 50 91 78 77 72 80 60 ...
##  $ Rape    : num  21.2 44.5 31 19.5 40.6 38.7 11.1 15.8 31.9 25.8 ...

prcomp(USArrests, scale = TRUE)

## Standard deviations (1, .., p=4):
## [1] 1.5748782744 0.9948694148 0.5971291155 0.4164493820
## 
## Rotation (n x k) = (4 x 4):
##                    PC1           PC2           PC3           PC4
## Murder   -0.5358994749  0.4181808654 -0.3412327280  0.6492278043
## Assault  -0.5831836349  0.1879856042 -0.2681484278 -0.7434074799
## UrbanPop -0.2781908746 -0.8728061931 -0.3780157931  0.1338777308
## Rape     -0.5434320914 -0.1673186354  0.8177779076  0.0890243227

prcomp(~ Murder + Assault + Rape, data = USArrests, scale = TRUE)

## Standard deviations (1, .., p=3):
## [1] 1.5357669768 0.6767948935 0.4282154425
## 
## Rotation (n x k) = (3 x 3):
##                   PC1           PC2           PC3
## Murder  -0.5826005597  0.5339532186 -0.6127565162
## Assault -0.6079818196  0.2140235835  0.7645600125
## Rape    -0.5393836250 -0.8179779129 -0.1999435899

plot(prcomp(USArrests))

summary(prcomp(USArrests, scale = TRUE))

## Importance of components%s:
##                             PC1       PC2       PC3       PC4
## Standard deviation     1.574878 0.9948694 0.5971291 0.4164494
## Proportion of Variance 0.620060 0.2474400 0.0891400 0.0433600
## Cumulative Proportion  0.620060 0.8675000 0.9566400 1.0000000

biplot(prcomp(USArrests, scale = TRUE))

Análisis factorial

El análisis factorial (AF) es tambien una técnica estadística de reducción de datos usada para explicar las correlaciones entre las variables observadas en términos de un número menor de variables no observadas llamadas factores.

El análisis factorial exploratorio, se utiliza para tratar de descubrir la estructura interna de un número relativamente grande de variables. La hipótesis a priori del investigador es que pueden existir una serie de factores asociados a grupos de variables. Las cargas de los distintos factores se utilizan para intuir la relación de éstos con las distintas variables. Es el tipo de análisis factorial más común.

El análisis factorial confirmatorio, AFC, trata de determinar si el número de factores obtenidos y sus cargas se corresponden con los que cabría esperar a la luz de una teoría previa acerca de los datos. La hipótesis a priori es que existen unos determinados factores preestablecidos y que cada uno de ellos está asociado con un determinado subconjunto de las variables. El análisis factorial confirmatorio entonces arroja un nivel de confianza para poder aceptar o rechazar dicha hipótesis.

El análisis de componentes principales es el método apropiado de extracción de factores, cuando el interés primordial se centra en la predicción o el número mínimo de factores necesarios para justificar la porción máxima de varianza representada en la serie de variables original, y cuando el conocimiento previo sugiere que la varianza específica y de error representan una porción relativamente pequeña de la varianza total. Por el contrario cuando se pretende identificar las dimensiones latentes o las construcciones representadas en las variables originales y se tiene poco conocimiento de la varianza específica y el error, lo más apropiado es utilizar el método factorial común. Si bien las complicaciones del análisis factorial común han contribuido al análisis generalizado de la técnica de componentes principales.

Para interpretar bien los factores se utiliza una rotación de ejes, ya que las soluciones factoriales no rotadas extraen factores según su orden de importancia. El primer factor tiende a ser un factor general por el que casi toda variable se ve afectada significativamente dando cuenta del mayor porcentaje de varianza. Los métodos de rotación ortogonales mas utilizados son VARIMAX, QUARTIMAX y EQUIMAX.

En la practica no se utiliza un criterio único de los factores a extraer, si bien se utiliza como primera aproximación el gráfico de autovalor para el criterio de contraste de caida, despues de interpretar los factores extraidos hay que valorar su caracter práctico.

Se leen la base de datos Seatbelts, que contiene datos mensuales de conductores de Gran Bretaña, muestos o heridos entre Enero de 1969 y Diciembre de 1983:

library(rela)
Belts <- Seatbelts[,1:7]
summary(Belts)

##  DriversKilled         drivers             front          
##  Min.   : 60.0000   Min.   :1057.000   Min.   : 426.0000  
##  1st Qu.:104.7500   1st Qu.:1461.750   1st Qu.: 715.5000  
##  Median :118.5000   Median :1631.000   Median : 828.5000  
##  Mean   :122.8021   Mean   :1670.307   Mean   : 837.2188  
##  3rd Qu.:138.0000   3rd Qu.:1850.750   3rd Qu.: 950.7500  
##  Max.   :198.0000   Max.   :2654.000   Max.   :1299.0000  
##       rear               kms           PetrolPrice        
##  Min.   :224.0000   Min.   : 7685.0   Min.   :0.08117889  
##  1st Qu.:344.7500   1st Qu.:12685.0   1st Qu.:0.09257748  
##  Median :401.5000   Median :14987.0   Median :0.10447669  
##  Mean   :401.2083   Mean   :14993.6   Mean   :0.10362400  
##  3rd Qu.:456.2500   3rd Qu.:17202.5   3rd Qu.:0.11405551  
##  Max.   :646.0000   Max.   :21626.0   Max.   :0.13302742  
##    VanKilled        
##  Min.   : 2.000000  
##  1st Qu.: 6.000000  
##  Median : 8.000000  
##  Mean   : 9.057292  
##  3rd Qu.:12.000000  
##  Max.   :17.000000

Calculo la matriz de correlaciones:

correl=cor(Belts,use="pairwise.complete.obs") 
correl

##               DriversKilled       drivers         front          rear
## DriversKilled  1.0000000000  0.8888263684  0.7067596370  0.3533510185
## drivers        0.8888263684  1.0000000000  0.8084114058  0.3436684967
## front          0.7067596370  0.8084114058  1.0000000000  0.6202247611
## rear           0.3533510185  0.3436684967  0.6202247611  1.0000000000
## kms           -0.3211015934 -0.4447631444 -0.3573823377  0.3330068886
## PetrolPrice   -0.3866060898 -0.4576675397 -0.5392393740 -0.1326272138
## VanKilled      0.4070411620  0.4853994918  0.4724207231  0.1217580839
##                         kms   PetrolPrice     VanKilled
## DriversKilled -0.3211015934 -0.3866060898  0.4070411620
## drivers       -0.4447631444 -0.4576675397  0.4853994918
## front         -0.3573823377 -0.5392393740  0.4724207231
## rear           0.3330068886 -0.1326272138  0.1217580839
## kms            1.0000000000  0.3839003754 -0.4980355895
## PetrolPrice    0.3839003754  1.0000000000 -0.2885584074
## VanKilled     -0.4980355895 -0.2885584074  1.0000000000

La prueba de esfericidad de Bartlett contrasta la hipótesis nula de que la matriz de correlaciones es una matriz identidad, en cuyo caso no existirían correlaciones significativas entre las variables y el modelo factorial no sería pertinente.

library(psych)
cortest.bartlett(Belts)

## R was not square, finding R from data

## $chisq
## [1] 968.7578864
## 
## $p.value
## [1] 1.260100725e-191
## 
## $df
## [1] 21

Para calcular los componentes principales basados en la matriz de correlaciones:

## Importance of components:
##                             Comp.1       Comp.2       Comp.3        Comp.4
## Standard deviation     1.920351456 1.2141728520 0.8489512909 0.78140438282
## Proportion of Variance 0.526821388 0.2106022449 0.1029597563 0.08722754421
## Cumulative Proportion  0.526821388 0.7374236329 0.8403833893 0.92761093348
##                               Comp.5        Comp.6         Comp.7
## Standard deviation     0.57544454265 0.32988394429 0.258386584886
## Proportion of Variance 0.04730520309 0.01554620239 0.009537661036
## Cumulative Proportion  0.97491613658 0.99046233896 1.000000000000
## 
## Loadings:
##               Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7
## DriversKilled -0.444         0.136  0.521 -0.360 -0.519  0.323
## drivers       -0.480         0.100  0.377         0.505 -0.596
## front         -0.475 -0.198        -0.103  0.417  0.414  0.612
## rear          -0.226 -0.686        -0.322  0.277 -0.406 -0.361
## kms            0.276 -0.627               -0.606  0.365  0.168
## PetrolPrice    0.327 -0.141  0.829  0.277  0.311              
## VanKilled     -0.334  0.259  0.526 -0.627 -0.387

A continuación se realiza un Analisis factorial con la libreria : psych, representando el gráfico de autovalor para el criterio de contraste de caida.

library(psych)
fa.parallel(Belts)

## The estimated weights for the factor scores are probably incorrect.  Try a different factor extraction method.
## The estimated weights for the factor scores are probably incorrect.  Try a different factor extraction method.

## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate =
## rotate, : A loading greater than abs(1) was detected. Examine the loadings
## carefully.

## The estimated weights for the factor scores are probably incorrect.  Try a different factor extraction method.

## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate =
## rotate, : An ultra-Heywood case was detected. Examine the results carefully

## The estimated weights for the factor scores are probably incorrect.  Try a different factor extraction method.

## Parallel analysis suggests that the number of factors =  3  and the number of components =  2

Belts.ps <- fa(Belts, fm="minres", nfactors=3, rotate="varimax")

## The estimated weights for the factor scores are probably incorrect.  Try a different factor extraction method.

## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate =
## rotate, : An ultra-Heywood case was detected. Examine the results carefully

print(Belts.ps) # Factor loadings as $loadings

## Factor Analysis using method =  minres
## Call: fa(r = Belts, nfactors = 3, rotate = "varimax", fm = "minres")
## Standardized loadings (pattern matrix) based upon correlation matrix
##                 MR1   MR3   MR2   h2       u2 com
## DriversKilled  0.81  0.27  0.26 0.80  0.19785 1.4
## drivers        0.88  0.39  0.25 1.00  0.00454 1.6
## front          0.47  0.54  0.71 1.01 -0.01177 2.7
## rear           0.21 -0.12  0.82 0.72  0.27655 1.2
## kms           -0.19 -0.93  0.32 1.00  0.00012 1.3
## PetrolPrice   -0.27 -0.44 -0.20 0.30  0.69660 2.1
## VanKilled      0.30  0.49  0.11 0.34  0.65652 1.8
## 
##                        MR1  MR3  MR2
## SS loadings           1.89 1.83 1.45
## Proportion Var        0.27 0.26 0.21
## Cumulative Var        0.27 0.53 0.74
## Proportion Explained  0.37 0.35 0.28
## Cumulative Proportion 0.37 0.72 1.00
## 
## Mean item complexity =  1.7
## Test of the hypothesis that 3 factors are sufficient.
## 
## The degrees of freedom for the null model are  21  and the objective function was  5.16 with Chi Square of  968.76
## The degrees of freedom for the model are 3  and the objective function was  0.07 
## 
## The root mean square of the residuals (RMSR) is  0.02 
## The df corrected root mean square of the residuals is  0.04 
## 
## The harmonic number of observations is  192 with the empirical chi square  1.85  with prob <  0.6 
## The total number of observations was  192  with Likelihood Chi Square =  13.41  with prob <  0.0038 
## 
## Tucker Lewis Index of factoring reliability =  0.922
## RMSEA index =  0.137  and the 90 % confidence intervals are  0.067 0.212
## BIC =  -2.36
## Fit based upon off diagonal values = 1

Analsis factorial con factanal

Belts.fa<-factanal(Belts,factors=3, rotation="varimax",scores = "Bartlett")
Belts.fa

## 
## Call:
## factanal(x = Belts, factors = 3, scores = "Bartlett", rotation = "varimax")
## 
## Uniquenesses:
## DriversKilled       drivers         front          rear           kms 
##         0.197         0.005         0.005         0.227         0.051 
##   PetrolPrice     VanKilled 
##         0.660         0.645 
## 
## Loadings:
##               Factor1 Factor2 Factor3
## DriversKilled  0.820   0.264   0.244 
## drivers        0.883   0.391   0.249 
## front          0.485   0.531   0.691 
## rear           0.204  -0.121   0.846 
## kms           -0.191  -0.903   0.311 
## PetrolPrice   -0.241  -0.469  -0.249 
## VanKilled      0.301   0.507         
## 
##                Factor1 Factor2 Factor3
## SS loadings      1.915   1.813   1.482
## Proportion Var   0.274   0.259   0.212
## Cumulative Var   0.274   0.533   0.744
## 
## Test of the hypothesis that 3 factors are sufficient.
## The chi square statistic is 11.54 on 3 degrees of freedom.
## The p-value is 0.00914

Análisis factorial con libreria : rela

library(rela)
res <- paf(as.matrix(Belts))
summary(res) # Automatically calculate KMO with MSA, determine the number of factors,

## $KMO
## [1] 0.66899
## 
## $MSA
##                   MSA
## DriversKilled 0.75338
## drivers       0.72295
## front         0.65795
## rear          0.40214
## kms           0.53517
## PetrolPrice   0.83413
## VanKilled     0.88197
## 
## $Bartlett
## [1] 968.76
## 
## $Communalities
##               Initial Communalities Final Extraction
## DriversKilled               0.79824          0.66181
## drivers                     0.87486          0.85141
## front                       0.87225          0.91049
## rear                        0.76216          0.69016
## kms                         0.67184          0.89761
## PetrolPrice                 0.36120          0.29405
## VanKilled                   0.36294          0.34590
## 
## $Factor.Loadings
##                   [,1]      [,2]
## DriversKilled  0.80934 -0.082307
## drivers        0.92263 -0.012744
## front          0.92638 -0.228709
## rear           0.41852 -0.717642
## kms           -0.53335 -0.783039
## PetrolPrice   -0.53096 -0.110142
## VanKilled      0.55008  0.208105
## 
## $RMS
## [1] 0.050492

# calculate chi-square of Bartlett's sphericity test, communalities and
# factor loadings. Communalities are 1 minus uniquenesses.
barplot(res$Eigenvalues[,1]) # First column of eigenvalues.

resv <- varimax(res$Factor.Loadings) # Varimax rotation is possible later.
print(resv)

## $loadings
## 
## Loadings:
##               [,1]   [,2]  
## DriversKilled  0.773 -0.254
## drivers        0.898 -0.211
## front          0.856 -0.422
## rear           0.255 -0.791
## kms           -0.689 -0.650
## PetrolPrice   -0.542       
## VanKilled      0.582       
## 
##                 [,1]  [,2]
## SS loadings    3.309 1.343
## Proportion Var 0.473 0.192
## Cumulative Var 0.473 0.664
## 
## $rotmat
##         [,1]     [,2]
## [1,] 0.97667 -0.21474
## [2,] 0.21474  0.97667

barplot(sort(colSums(loadings(resv)^2),decreasing=TRUE)) # screeplot using rotated SS loadings.

scores <- as.matrix(Belts) %*% as.matrix(resv$loadings) # Get factor scores in a simple manner.

El índice KMO compara los valores de las correlaciones entre las variables y sus correlaciones parciales. Si el índice KMO está próximo a 1, el ACP se puede hacer. Si el índice es bajo (próximo a 0), el ACP no será relevante.

Análisis cluster

Introducción

El análisis cluster (AC) es un conjunto de técnicas multivariantes cuyo principal propósito es agrupar objetos basándose en las características que poseen. El AC clasifica los objetos en clases ó conglomerados de tal forma que cada objeto sea parecido a los que hay en el conjunto de su conglomerado. Los conglomerados resultantes tendrán que tener un alto grado de homogeneidad interna (dentro del conglomerado) y de heterogeneidad externa (entre conglomerados).

Para realizar un AC se necesita una medida de similitud, de correspondencia, o parecido entre objetos que van a ser analizados. Esta puede ser medida de varias formas, pero en el AC, tres son los métodos de simulitud que se utilizan: medidas de correlación (la matriz de correlación entre variables), medidas de distancia (la más utilizada es la distancia euclidea) y las medidas de asociación, cuando se tienen datos no métricos.

Una vez se tiene una matriz de simulitud calculada, hay que realizar el proceso de partición de los datos, para ello hay que decidir el algoritmo de aglomeración utilizado en la formación de conglomerados, para despues tomar una decisión acerca del número de conglomerados que se van a formar. Los algoritmos de obtención de conglomerados se subdividen en métodos jerárquicos (encadenamiento simple, encadenamiento completo, encadenamiento médio, método de Ward, método del centroide), de los no jerarquicos (umbral secuencial, umbral paralelo, optimización), o una combinación de los dos.

El dendrograma ó gráfico en forma de arbol, es una herramienta visual que ayudar a decidir el número de conglomerados que podrían representar mejor la estructura de los datos.

Para determinar el número final de conglomerados a formar o regla de parada, no hay un procedimiento determinado, ha de decidirlo el investigador en la fase de interpretación de los datos.

Cluster Jerarquico

R incluye la función hclust para la realización de clasificaciones jerárquicas.

Una clasificación precisa que el investigador le pase dos opciones básicas a R:

●Medida de similitud/disimilitud para calcular la matriz de distancias

●Estrategia de fusión

Posteriormente, en la interpretación de los datos el investigador deberá decidir el punto de corte y, con ello, el número de grupos.

Se puede usar la función nativa de R dist para calcular la matriz de distancias; por defecto el programa usa la distancia euclídea (euclidean), siendo posible elegir varias más con la opción method=“”: maximum, manhattan, canberra, binary, minkowski.

Las estrategias de fusión disponibles en hclust son las siguientes:

●“ward”: la distancia entre dos conglomerados es la suma de los cuadrados entre dos conglomerados sumados para todas las variables. En cada paso del procedimiento de aglomeración se minimiza la suma de cuadrados dentro del conglomerado para todas las particiones. Tiende a combinar los conglomerados con un número reducido de dimensiones.

●“single”: distancia mínima o vecino más próximo, encuentra dos objetos separados por la distancia más corta y los coloca en un primer conglomerado. A continuación se encuetra la distancia más corta, y se une a el un tercer objeto o se forma un nuevo conglomerado de dos miembros. El proceso continua hasta que todos los objetos están en un conglomerado.

●“complete”: Parecido al “single”(simple) pero utilizando como criterio de agregación la distancia máxima (aproximación al vecino más lejano).

●“average”: comienza igual que los métodos “single” y “complete”, pero el método de aproximación es la distancia media entre todos los individuos de un conglomerado y los de otro.

●“mcquitty”

●“median”

●“centroid”: la distancia entre dos conglomerados es la distancia (euclidia) entre centroides.

Realizamos un Cluster Jerarquico con las tasas de crimen en USA (USArrests), utilizando una medida de similaridad de distancias euclideas y el método de agrupación de encadenamiento medio:

require(graphics)
str(USArrests)

## 'data.frame':    50 obs. of  4 variables:
##  $ Murder  : num  13.2 10 8.1 8.8 9 7.9 3.3 5.9 15.4 17.4 ...
##  $ Assault : int  236 263 294 190 276 204 110 238 335 211 ...
##  $ UrbanPop: int  58 48 80 50 91 78 77 72 80 60 ...
##  $ Rape    : num  21.2 44.5 31 19.5 40.6 38.7 11.1 15.8 31.9 25.8 ...

hc <- hclust(dist(USArrests), "ave")
plot(hc)

plot(hc, hang = -1)

Se realiza el mismo cluster Jerarquico, se establece una parada de 10 cluster (cutree), y se realiza un nuevo cluster con los 10 conglomerados.

hc <- hclust(dist(USArrests)^2, "cen")
memb <- cutree(hc, k = 10)
cent <- NULL
for(k in 1:10){
  cent <- rbind(cent, colMeans(USArrests[memb == k, , drop = FALSE]))
}
str(cent)

##  num [1:10, 1:4] 11.47 13.5 9.95 11.5 5.59 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : chr [1:4] "Murder" "Assault" "UrbanPop" "Rape"

hc1 <- hclust(dist(cent)^2, method = "cen", members = table(memb))
opar <- par(mfrow = c(1, 2))
plot(hc,  labels = FALSE, hang = -1, main = "Original Tree")
plot(hc1, labels = FALSE, hang = -1, main = "Re-start from 10 clusters")

par(opar)

Determinar el número de cluster con k-means

Determinar el número de cluster o la regla de parada es una de las tareas más sensibles en el AC, una forma de incluir un test para decidir el numero de conglomerados en utilizar la función k-meas de R, y representar gráficamente el vector la suma de cuadrados entre conglomerados (withniss), para diferentes agrupamientos.

d=dist(USArrests)^2

# Determinar el numero de cluster

# Fija semilla inicial de numeros pseudoaleatorios
# para obtener misma serie aleatoria en ejecuciones.
set.seed(123)
 
# Crea vector "Errores", sin datos
# Crea variable "K_Max" con la cant. maxima de k a analizar
Errores <-NULL
K_Max   <-10
 

# Ejecuta kmeans con diferentes cluster, desde 1 hasta 10
# Luego guarda el error de cada ejecucion en el vector "Errores"
for (i in 1:K_Max)
  {
  Errores[i] <- sum(kmeans(d[-1], centers=i)$withinss)
  }
# Grafica el vector "Errores"
plot(1:K_Max, Errores, type="b", 
              xlab="Cantidad de Cluster", 
              ylab="Suma de error")

Cluster Jerárquico de series temporales con DTW

La agrupación en series temporales consiste, como se apreció en el apartardo anterior, en dividir el data-set de series temporales en grupos basados en similitudes o distancias, de modo que las series temporales de un mismo grupo sean similares.

Para la agrupación de series de tiempo, el primer paso es elaborar una métrica de distancia ó similitud apropiada, y luego, en el segundo paso, utilizar técnicas de agrupación existentes, como k-means, clustering jerárquico, clustering basado en densidad o clustering subespacial.

Para encontrar estructuras de agrupamiento, de series temporales una metrica adecuada es la distancia Dynamic Time Warping (DTW), que encuentra la alineación óptima entre dos series temporales.

El Alineamiento Temporal Dinámico (Dynamic Time Warping, DTW), es una técnica desarrollada en el ambito de reconocimiento de la voz, cuyo objetivo es comparar dos emisiones vocales aisladas a fin de determinar si corresponden o no a la misma palabra.

De manera resumida, el DTW se diseño para comprobar si existe una sincronización temporal (alineamiento temporal) entre dos grupos fónicos. La falta de alineamiento, por lo general, no obedece a una ley fija (p. e., un retardo constante), sino que se da de forma heterogénea, produciéndose así variaciones localizadas que aumentan o disminuyen la duración del tramo de análisis.

El problema de calcular una función de alineamiento se resuelve en el DTW mediante técnicas de programación dinámica, motivo por el que se conoce a esta técnica como alineamiento temporal dinámico.

En el ejemplo, se utiliza un conjunto de datos de la serie temporal de gráficos de control sintético, que contiene 600 ejemplos de gráficos de control. Cada gráfico de control es una serie de tiempo con 60 valores. Hay seis clases:

1. 1-100 Normal,

2. 101-200 Cíclico,

3. 201-300 Tendencia creciente,

4. 301-400 Tendencia decreciente,

5. 401-500 Cambio hacia arriba y

6. 501-600 Descenso cambio.

Los datos y el script se ha obtenido de :http://www.rdatamining.com/examples/time-series-clustering-classification.

sc <- read.table("http://kdd.ics.uci.edu/databases/synthetic_control/synthetic_control.data", header=F, sep="")
# randomly sampled n cases from each class, to make it easy for plotting
n <- 10
s <- sample(1:100, n)
idx <- c(s, 100+s, 200+s, 300+s, 400+s, 500+s)
sample2 <- sc[idx,]
observedLabels <- c(rep(1,n), rep(2,n), rep(3,n), rep(4,n), rep(5,n), rep(6,n))
# compute DTW distances
library(dtw)

## Loading required package: proxy

## 
## Attaching package: 'proxy'

## The following objects are masked from 'package:stats':
## 
##     as.dist, dist

## The following object is masked from 'package:base':
## 
##     as.matrix

## Loaded dtw v1.18-1. See ?dtw for help, citation("dtw") for use in publication.

distMatrix <- dist(sample2, method="DTW")
# hierarchical clustering
hc <- hclust(distMatrix, method="average")
plot(hc, labels=observedLabels, main="")

Cluster Jerárquico con span-tree

Si al clasificar un conjunto de datos no encontramos una decisión adecuada para establecer una regla de parada, y obtenemos como resultado un conjunto inclasificable, es util combinar la metodología de jerarquía y topología de los árboles de expansión mínima (MST, Minimal Spanning Tree).

ado un grafo, su árbol mínimo generador (o árbol de peso mínimo o árbol mínimo de expansión) es un árbol que pasa por todos los vértices y que la suma de sus aristas es la de menor peso.

La gráfica siguiente ilustra el proceso, se trata de un arbol de expansión de 7 nodos (acciones), con las medidas de distancia que hay entre ellos. Partiendo de la distancia menor, 29, entre el nodo 8 y 4, se inicia el arbol, que en el nodo 4, encuentra la segunda minima distancia en su enlace con el nodo 2, 31, y en el nodo 2, la menor distancia la encuentra con el nodo 5, operando de esta manera hasta tener enlazados todos los nodos. Como se puede observar en la figura, un nodo puede tener, más de un enlace con otros nodos, tal y como ocurre en el gráfico, nodo 3.

Curva ROC

Basandonos en el trabajo pionero de Mantegna (1999) sobre tasas de retornos en los mercados financieros, realizamos un ejercicio clasificatorio con la función span-tree, de la librería vegan, que ordena y establece dependencias entre diversos activos financieros utilizando datos de rendimientos mensuales de activos procedentes de Berndt’s The Practice of Econometrics: http://web.stanford.edu/~clint/berndt/

# read prices from csv file
bolsa.df = load("berndt.RData")
colnames(berndt.df)

##  [1] "Año"    "MOBIL"  "TEXACO" "IBM"    "DEC"    "DATGEN" "CONED" 
##  [8] "PSNH"   "WEYER"  "BOISE"  "MOTOR"  "TANDY"  "PANAN"  "DELTA" 
## [15] "CONTIL" "CITGRP" "GERBER" "GENMIL" "MARKET" "RKFREE"

#[1] "Año"    "MOBIL"  "TEXACO" "IBM"    "DEC"    "DATGEN" "CONED"  "PSNH"   "WEYER"  "BOISE" 
#[11] "MOTOR"  "TANDY"  "PANAN"  "DELTA"  "CONTIL" "CITGRP" "GERBER" "GENMIL" "MARKET" "RKFREE

# create zooreg object - regularly spaced zoo object


library(zoo)
berndt.z = zooreg(berndt.df[,-1], start=c(1978, 1), end=c(1987,12),frequency=12)
index(berndt.z) = as.yearmon(index(berndt.z))

start(berndt.z)

## [1] "ene. 1978"

end(berndt.z)

## [1] "dic. 1987"

nrow(berndt.z)

## [1] 120

# note: coredata() function extracts data from zoo object
returns.mat = as.matrix(coredata(berndt.z))

# create the correlation coefficients

coef.corr <- cor(returns.mat)
coef.d <- (1-coef.corr^2) # compute distance (Mantegna, 1998)

# hierarchical cluster whir hclust

d <- as.dist(as.matrix(coef.d)) # find distance matrix 

#Function spantree finds a minimum spanning tree (MTA) connecting all points, but disregarding dissimilarities that are at or above the threshold or NA.


library(vegan)

## Loading required package: permute

## Loading required package: lattice

## This is vegan 2.4-3

tr <- spantree(d)
## Add tree to a metric scaling 
plot(tr, cmdscale(d), type = "t")

## Find a configuration to display the tree neatly
plot(tr, type = "t")

## Initial stress        : 0.10907
## stress after   2 iters: 0.06448

## Depths of nodes
depths <- spandepth(tr)
plot(tr, type = "t", label = depths)

## Initial stress        : 0.10907
## stress after   2 iters: 0.06448

## Plot as a dendrogram
cl <- as.hclust(tr)
plot(cl)

Inferencia no parámetrica

Introducción

Se dice que se ajusta el modelo paramétrico cuando se estiman sus parámetros a partir de un conjunto de observaciones que siguen dicho modelo, de manera que pueden hacerse predicciones de nuevos valores de \(Y\) conocido el valor de \(X\), y tener información precisa acerca de la incertidumbre asociada a la estimación y a la predicción. Sin embargo, si el modelo paramétrico no es el adecuado al análisis de datos que estamos realizando, pueden llevar a conclusiones que queden muy alejadas de la realidad, dado que el modelo paramétrico conlleva un grado de exactitud en las afirmaciones que de él se derivan y que son adecuadas siempre y cuando se cumplan los supuestos básicos sobre los que se apoya su construcción teórica. De hecho, los modelos paramétricos presentan una estructura teórica tan rígida que no pueden adaptarse a muchos conjuntos de datos de los que hoy día se disponen para el análisis económico.

La econometría no paramétrica aparece como consecuencia de intentos por solucionar problemas que existen en la econometría paramétrica como, por ejemplo, la consistencia entre los datos y los principios de maximización, homocedasticidad, o la necesidad de asumir una determinada relación, por lo general de forma lineal entre las variables de interés. Esta preocupación llevó a una serie de investigadores a utilizar formas funcionales flexibles para aproximarse a relaciones desconocidas entre las variables. El plantear formas funcionales flexibles requiere el conocimiento del valor esperado de la variable \(Y\), condicional en las otras, \(X\). Esto conlleva la necesidad de estimar la función de densidad de \(Y\) condicional en \(X\).

La econometría no paramétrica no parte de supuestos sobre la distribución de probabilidad de las variables bajo estudio, sino que trata de estimar dicha distribución para encontrar la media condicional y los momentos de orden superior (por ejemplo, la varianza) de la variable de interés. Una de las desventajas de este método es, sin embargo, la necesidad de contar con muestras muy grandes si es que se desea estimar la función de relación entre ambas variables de manera precisa. Además el tamaño de la muestra debe aumentar considerablemente conforme aumenta el número de variables involucradas en la relación.

Los modelos de regresión paramétricos suponen que los datos observados provienen de variables aleatorias cuya distribución es conocida, salvo por la presencia de algunos parámetros cuyo valor se desconoce.

\[\begin{equation} y_i= \beta_1 + \beta_2 x_i + e_i, i=1,2,…, n \end{equation}\]

.

con \(e_i\approx N(0;\sigma^2)\)

Este es un modelo estadístico con tres parámetros desconocidos:\(\beta_1\); \(\beta_2\) y \(\sigma^2\).

Una formulación general de un modelo de regresión paramétrico es la siguiente:

\(y_i=m(x_i;\theta)+e_i, i=1,2,…, n\), \(\theta \in \Theta \in \Re\)

Donde \(m(x_i;\theta)\) es una función conocida de \(X\) y de \(\theta\), que es desconocido, \(e_i\) es una variable aleatoria idénticamente distribuida con \(E(e_i)=0\) y \(Var(e_i)=\sigma^2\). El modelo de regresión lineal simple sería un caso particular con \(\theta=(\beta_1, \beta_2)\) y \(m(x_i;\theta)=\beta_1 + \beta_2 x_i)\).

Donde los valores de la variable explicativa son conocidos, por lo que se dice que el modelo tiene diseño fijo, y dado que la varianza de los errores es constante el modelo es Homocedástico.

Considerando una variable aleatoria bivariante con densidad conjunta \((Y,X)\) , cabe definir la función de regresión \(f(y,x)\) como \(m(x,\theta)=E(Y/X=x)\), es decir el valor esperado de \(Y\) cuando \(X\) toma el valor conocido \(x\) . Entonces \(E(Y/X)=m(X,\theta))\), y definiendo \(e=Y-m(X,\theta)\), se tiene que: \(Y=m(X,\theta)+e\) ,\(E(e/X)=0\) y \(Var(e/X)=\sigma^2\)

Se supone que se observan ahora \(n\) pares de datos \((y_i,x_i),i=1…n\). Estos datos siguen el modelo de regresión no paramétrico:\(y_i=m(x_i;\theta)+e_i\).

Una vez establecido el modelo, el paso siguiente consiste en estimarlo (o ajustarlo) a partir de las n observaciones disponibles. Es decir hay que construir un estimador de la función de regresión y un estimador de la varianza del error. Los procedimientos de estimación de se conocen como métodos de suavizado.

El abanico de técnicas disponibles para estimar no paramétricamente la función de regresión es amplísimo e incluye, entre otras, las siguientes:

• Ajuste local de modelos paramétricos. Se basa en hacer varios (o incluso infinitos, desde un punto de vista teórico) ajustes paramétricos teniendo en cuenta únicamente los datos cercanos al punto donde se desea estimar la función.

• Suavizado mediante splines. Se plantea el problema de buscar la función \(\hat m(x)\) que minimiza la suma de los cuadrados de los errores (\(e_i=y_i-\hat m(x)\) ) más un término que penaliza la falta de suavidad de las funciones (\(\hat m(x)\) ) candidatas (en términos de la integral del cuadrado de su derivada segunda).

• Métodos basados en series ortogonales de funciones. Se elige una base ortonormal del espacio vectorial de funciones y se estiman los coeficientes del desarrollo en esa base de la función de regresión. Los ajustes por series de Fourier y mediante wavelets son los dos enfoques más utilizados.

• Técnicas de aprendizaje supervisado. Las redes neuronales, los k vecinos más cercanos y los árboles de regresión se usan habitualmente para estimar \(m(x)\).

Ajuste local (estimadores nucleo).

Los histogramas son siempre, por naturaleza, funciones discontinuas; sin embargo, en muchos casos es razonable suponer que la función de densidad de la variable que se está estimando es continua. En este sentido, los histogramas son estimadores insatisfactorios. Los histogramas tampoco son adecuados para estimar las modas, a lo sumo, pueden proporcionar “intervalos modales”, y al ser funciones constantes a trozos, su primera derivada es cero en casi todo punto, lo que les hace completamente inadecuados para estimar la derivada de la función de densidad.

Los estimadores de tipo núcleo (o kernel) fueron diseñados para superar estas dificultades. La idea original es bastante antigua y se remonta a los trabajos de Rosenblatt y Parzen en los años 50 y primeros 60. Los estimadores kernel son, sin duda, los más utilizados y mejor estudiados en la teoría no paramétrica.

Dada una m.a.s (\(X_1,X_2,...X_n\)) con densidad \(f\) , estimamos dicha densidad en un punto \(t\) por medio del estimador:

\[\hat f(t)=\frac{1}{nh}\sum_{i=1}^n K(\frac{t-X_i}{h})\]

donde \(h\) es una sucesión de parámetros de suavizado, llamados ventanas o amplitudes de banda (windows, bandwidths) que deben tender a cero ”lentamente" (\(h \Rightarrow 0\) ,\(nh\Rightarrow inf\) ) para poder asegurar que \(\hat f\) tiende a la verdadera densidad \(f\) de las variables \(X_i\) y \(K\) es una función que cumple \(\int K=1\). Por ejemplo:

• Nucleo Gaussiano \(\frac{1}{\sqrt{2 \pi}}e^{- \frac {u^2} 2}\)

• Núcleo Epanechnikov \(\frac 3 4 (1-u^2)I_{\mid u \mid < 1}\)

donde \(I_{\mid u \mid < 1}\) es la función que vale 1 si \({\mid u \mid < 1}\) y cero si \({\mid u \mid \geq 1}\)

• Núcleo Triangular \((1-\mid u \mid)I_{\mid u \mid < 1}\)

• Núcleo Uniforme \(\frac 1 2 I_{\mid u \mid < 1}\)

• Núcleo Biweight \(\frac {15} {16} (1-u^2)I_{\mid u \mid < 1}\)

• Núcleo Triweight \(\frac {35} {32} (1-u^2)I_{\mid u \mid < 1}\)

Los programas informáticos eligen el ancho de ventana, \(h\) siguiendo criterios de optimización (error cuadrático medio).

En R la estimación de una función de densidad kernel se realiza con la función “density”, con los datos del vector x hay que realizar el siguiente programa:

x <- c(2.1,2.6,1.9,4.5,0.7,4.6,5.4,2.9,5.4,0.2)
density(x,kernel="epanechnikov")

## 
## Call:
##  density.default(x = x, kernel = "epanechnikov")
## 
## Data: x (10 obs.);   Bandwidth 'bw' = 1.065
## 
##        x                 y         
##  Min.   :-2.9942   Min.   :0.0000  
##  1st Qu.:-0.0971   1st Qu.:0.0237  
##  Median : 2.8000   Median :0.0943  
##  Mean   : 2.8000   Mean   :0.0862  
##  3rd Qu.: 5.6971   3rd Qu.:0.1524  
##  Max.   : 8.5942   Max.   :0.1695

plot(density(x,kernel="epanechnikov"))

Los estimadores núcleo establecen que el peso de (\(X_i\),\(Y_i\)) en la estimación de \(m(x)\) es:

\[W_i(t,X_i)=\frac {\frac 1 h K(\frac{t-X_i}{h})}{\hat f(t)}\]

donde \(K(t)\) es una función de densidad simétrica (por ejemplo, la normal estándar) y \(\hat f(t)\) es un estimador kernel de la densidad como el definido en el apartado anterior.

\(W_i(t,X_i)\) es, para cada i, una función de ponderación que da “mayor importancia” a los valores \(X_i\) de la variable auxiliar que están cercanos a t.

Una expresión alternativa para \(W_i(t,X_i)\)

\[W_i(t,X_i)=\frac {K(\frac{t-X_i}{h})}{\sum_{j=1}^nK(\frac{t-X_i}{h})}\]

A partir de los pesos puede resolverse el problema de mínimos cuadrados ponderados siguiente:

\[min_{a,b}\sum_{i=1}^n W_i(Y_i-(\beta_0+\beta_1(t-X_i)))^2\]

los parámetros así obtenidos dependen de \(t\), porque los pesos \(W_i\) también dependen de \(t\), la recta de regresión localmente ajustada alrededor de \(t\) sería:

\(l_t(X)=\hat \beta_0(t)+\hat \beta_1(t)(t-X_i)\)

Y la estimación de la función en el punto en donde

\(\hat m(t) =l_t(t)=\hat \beta_0(t)\)

Las funciones núcleo usadas en la estimación no paramétrica de la regresión son las mismas que en la densidad.

Si se generaliza al ajuste local de regresiones polinómicas de mayor grado, es decir si pretendemos estimar una forma lineal del tipo:

\(\beta_0+\beta_1 X_i+\beta_3 X_i^2+...+\beta_q X_i^q\)

con la salvedad de que en vez del valor \(X_i\) en la regresión lineal múltiple se utiliza el valor \(t-X_i\) . El estimador de polinomios locales de grado \(q\) asignado los pesos obtenidos mediante la función núcleo se resuelve el siguiente problema de regresión polinómica ponderada:

\[min_{a,b}\sum_{i=1}^n W_i(Y_i-(\beta_0+\beta_1(t-X_i)+\beta_2(t-X_i)^2+...+\beta_q(t-X_i)^q))^2\]

Los parámetros \(\hat \beta_j= \hat \beta_j(t)\) dependen del punto t en donde se realiza la estimación, y el polinomio ajustado localmente alrededor de t sería:

\(P_{q,t}=\sum_{j=0}^q \hat \beta_j(t-X_i)^j\)

Siendo \(m(t)\) el valor de dicho polinomio estimado en el punto en donde \(X=t\):

\(\hat m(t) =P_{q,o}= \hat \beta_0(t)\)

En el caso particular del ajuste de un polinomio de grado cero, se obtiene el estimador de Nadaraya −Watson, o estimador núcleo de la regresión:

\[\hat m_k(t)=\frac {\sum_{i=0}^n K(\frac{t-X_i}{h})Y_i}{\sum_{i=0}^n K(\frac{t-X_i}{h})}\]

El estimador del parámetro de suavizado \(h\) tiene una importancia crucial en el aspecto y propiedades del estimador de función de regresión. Valores pequeños de dan mayor flexibilidad al estimador y le permiten acercarse a todos los datos observados, pero originan altos errores de predicción (sobre-estimación), valores mas altos de h ofrecerán un menor grado de ajustes a los datos pero predicican mejor, pero si h es demasiado elevado tendremos una falta de ajuste a los datos (sub-estimación).

Utilizando la base de datos “cars” de R, que contine las variables “dist” (distancia de parada) y “speed” (velocidad), vamos a realizar la representación gráfica de la regresión kernel realizada con el estimador de Nadaraya–Watson con diferentes parámetros de suavizado.

data(cars)
plot(cars$speed, cars$dist)
lines(ksmooth(cars$speed, cars$dist, "normal", bandwidth = 2), col = 2)
lines(ksmooth(cars$speed, cars$dist, "normal", bandwidth = 5), col = 3)

Splines

Para poder estimar la función \(f\) de la forma más sencilla posible, deberíamos poder representarla de forma que \(Y_i=f(x_i)+e_i\),\(i=1,2,...,n\) se convierta en un modelo lineal.

Y esto se puede hacer eligiendo una base de funciones de dimensión \(q\) que genere un subespacio de funciones que incluya a \(f\) como elemento y que pueda expresarse como:

\[f(x)=\sum_{j=1}^q \beta_j s_j(x)\]

Siendo \(\beta_j\) un parámetro desconocido, asociado al elemento \(j\), \(s_j(x)\) de dicha base de funciones.

De manera que:

\[Y_i=\sum_{j=1}^q \beta_j s_j(x) + e_i\] (1)

Se convierte en un modelo lineal de dimensión \(q\).

La regresión con funciones base polinómicas es la propuesta más sencilla para este tipo de estimaciones.

Supongamos que \(f\) es un polinomio de grado 4 de forma que el espacio de polinomios de grado 4 contiene a \(f\). Una base de este subespacio es:

\[\begin{equation} \left\lbrace\begin{array}{ll} s_1(x)=1\\ s_2(x)=x\\ s_3(x)=x^2\\ s_4(x)=x^3\\ s_5(x)=x^4 \end{array}\right.\end{equation}\]

De manera que el modelo (1) se convierte en:

\[Y_i=\beta_1 + \beta_2 x_i + \beta_3 x_i^2 + \beta_4 x_i^3 + \beta_5 x_i^4 + e_i\]

Un spline es una curva diferenciable definida en porciones mediante polinomios, que se utiliza como bases de funciones para aproximar curvas con formas complicadas.

Las bases de spilines más populares:

• Bases de polinomios truncados.

• Bases de splines cúbicos.

• Bases de B-splines.

• Bases de thin plate splines.

Una función spline está formada por varios polinomios, cada uno definido sobre un subintervalo, que se unen entre sí obedeciendo a ciertas condiciones de continuidad.

Supongamos que se ha fijado un entero \(q \neq 0\), de manera que disponemos de \(q+1\) puntos, a los que denominaremos nodos, tales que \(t_0<t_1<t_2<...<t_q\), en los que troceamos nuestro conjunto de datos. Decimos entonces que una función spline de grado \(q\) con nodos en \(t_0,t_1,t_2,...,t_q\) es una función \(S\) que satisface las condiciones:

1. en cada intervalo \([t_{j-1},t_j)\) , \(S\) es un polinomio de grado menor o igual a \(q\).

2. \(S\) tiene una derivada de orden \((q-1)\) continua en \([t_0,t_q ]\).

Los splines de grado \(0\) son funciones constantes por zonas. La expresión matemática de un spline de grado \(0\) es la siguiente:

\[\begin{equation} \left\lbrace\begin{array}{ll} s_0(x)=c_0 & x \in [t_0,t_1) \\ s_j(x)=c_j & x \in [t_{j-1},t_j) \\ ...\\ s_{q-1}(x)=c_{q-1} & x \in [t_{q-1},t_q) \end{array}\right.\end{equation}\]

En la figura siguiente se muestran las gráficas correspondientes a los splines de grado cero.

Splines de grado 0

Los splines de grado \(0\), se define en un solo tramo de nudo y ni siquiera es continua en los nudos. Equivale a realizar una regresión por tramos.

\[Y_i=\beta_0 c_0 x_i + \beta_1 c_1 x_i + ... + \beta_{q-1} c_{q-1} x_i+ e_i\]

donde \(c_j\) toma valor \(1\) en \([t_{j-1},t_j)\) y cero en el resto.

Un spline de grado 1 o lineal se expresa como:

\[\begin{equation} \left\lbrace\begin{array}{ll} s_0(x)=a_0x + b_0 & x \in [t_0,t_1) \\ s_j(x)=a_jx + b_j & x \in [t_{j-1},t_j) \\ ...\\ s_{q-1}(x)=a_{q-1}x + b_{q-1} & x \in [t_{q-1},t_q) \end{array}\right.\end{equation}\]

La representación gráfica de un spline lineal aparece en la figura siguiente:

Splines lineal

Las funciones de spilines más comúnmente utilizadas son las de grado 3 ó cúbicas. Son polinomios de grado tres a trozos, que son continuos en los nodos al igual que su primera y segunda derivada, proporcionando un excelente ajuste a los puntos tabulados y a través de cálculo que no es excesivamente complejo.

Sobre cada intervalo \([t_{j-1},t_j)\), \(S\) está definido por un polinomio cúbico diferente.

\[\begin{equation} \left\lbrace\begin{array}{ll} s_0(x)=a_0x^3 + b_0x^2 + c_0x + d_0 & x \in [t_0,t_1) \\ s_j(x)=a_jx^3 + b_jx^2 + c_jx + d_0 & x \in [t_{j-1},t_j) \\ ...\\ s_{q-1}(x)=a_{q-1}x^3 + b_{q-1}x^2 + c_{q-1}x + d_{q-1} & x \in [t_{q-1},t_q) \end{array}\right.\end{equation}\]

Los polinomios\(S_{j-1}\) y \(S_j\) interpolan el mismo valor en el punto \(t_j\), para que se cumpla:

\(S_{j-1}(x_i)=y_i=S_j(x_i)\)

por lo que se garantiza que \(S\) es continuo en todo el intervalo. Además, se supone que \(S'\) y \(S''\) son continuas, condición que se emplea en la deducción de una expresión para la función del spline cúbico.

Aplicando las condiciones de continuidad del spline \(S\) y de las derivadas primera \(S\) y segunda \(S''\), es posible encontrar la expresión analítica del spline.

Una de las bases de splines cúbicos más utilizadas basadas en \(q-2\) nodos interiores, \(x^*_j\), \(j=1,2,...,q-2\) es:

\[\begin{equation} \left\lbrace\begin{array}{ll} s_0(x)=1 \\ s_1(x)=x \\ s_{j+2}(x)=R[x_j,x^*_j] \end{array}\right.\end{equation}\]

Siendo \(R[x,z]=\frac 1 4[(z-\frac 1 2)^2-\frac 1 {12}][(x-\frac 1 2)^2-\frac 1 {12}]-\frac 1 {24}[(\mid z- x \mid - \frac 1 2)^4- [(\mid z- x \mid - \frac 1 2)^2+\frac 7 {240}]\)

Con esta base de splines definimos \(f\) a través de un modelo lineal con matriz de regresores \(X\) con \(n\) filas y \(q\) columnas cuya i_esima fila es:

\(X_i=[1,x_i,R[x_i,x^*_1],R[x_i,x^*_2],...,R[x_i,x^*_{q-2}]]\)

Los elementos de una base de splines cúbicos son polinomios de grado 3. Un Spline cúbico se representa en la figura siguiente:

Splines lineal

Un tema importante es la elección del grado de suavización del spline. Una de las posibilidades es a través del contraste de hipótesis, valorar la posibilidad de utilizar uno o más nodos. Pero lo aconsejado es mantener fija la base de splines y controlar el grado de suavización añadiendo una penalización a la función objetivo de mínimos cuadrados:

\(\lambda \beta' S \beta\)

Donde \(S\) es una matriz de orden \(q*q\) con coeficientes conocidos que dependen de la base elegida y un parámetro de suavizado \(\lambda\).

La solución del modelo de regresión lineal penalizado en donde la matriz de regresores está ahora definida por la base de splines y la penalización sería:

\(\hat \beta = (X'X - \lambda S)^{-1}X'y\)

El modelo de regresión lineal con spilines penalizados es equivalente al siguiente modelo de regresión lineal:

\(Y=\beta X + e\)

En donde \(Y'=(y,0,0...0)\) es un vector de dimensión \((n+q)\), es decir el vector \(y\) seguido de tantos ceros como nodos se han utilizado en la base de los splines.

La matriz de regresores

\[X' = \begin{bmatrix}X\\ B \sqrt \lambda\\ \end{bmatrix}\]

tiene ahora orden \((n+q)*q\), siendo \(B\) una matriz que cumple \(B=S'S\) y que se obtiene a través de la descomposición de Cholesky, \(\lambda\) el parámetro de suavizado y \(e\) un vector de \((n+q)\) errores aleatorios.

El parámetro de suavización,\(\lambda\),es a priori desconocido y hay que determinarlo, si es muy alto suaviza los datos en exceso, un criterio utilizado para elegir el parámetro es del valor que minimiza el estadístico general de validación cruzada.

La regresión por splines puede realizarse con múltiples variables explicativas, si tenemos ahora dos explicativas, \(x_i\) y \(z_i\),y queremos estimar el siguiente modelo aditivo:

\(y_i=f_1(x_i)+f_2(z_i)+e_i\)

Representaríamos cada una de estas dos funciones a través de una base de splines penalizados, que tomando la base cúbica quedaría:

\(f_1(x)=\delta_1+\delta_2X_i+\sum_{j=1}^{q-2} \delta_jR[x_i,x_j^*]\)

\(f_2(z)=\gamma_1+\gamma_2X_i+\sum_{j=1}^{q-2} \gamma_jR[z_i,z_j^*]\)

Partiendo de la base de datos “cars”, la función R “smooth.spline” realiza la regresión por splines utilizando una base de splinee cúbicos penalizados:

attach(cars)
plot(speed, dist, main = "data(cars) &  smoothing splines")
cars.spl1 <- smooth.spline(speed, dist)
cars.spl1

## Call:
## smooth.spline(x = speed, y = dist)
## 
## Smoothing Parameter  spar= 0.78013  lambda= 0.11122 (11 iterations)
## Equivalent Degrees of Freedom (Df): 2.6353
## Penalized Criterion (RSS): 4187.8
## GCV: 244.1

cars.spl2 <- smooth.spline(speed, dist,spar=0.10)
lines(cars.spl1, col = "blue")
lines(cars.spl2, col = "red")

En la función “smooth.spline” el parámetro de suavizado es un valor generalmente entre 0 y 1, en tanto que el coeficiente \(\lambda\) se obtiene en el criterio de aceptación (logaritmo de verosimilitud penalizado). En el ejercicio el programa elige un \(\lambda=0.7801305\). Si se desea un función menos suavizada habrá que elegir un parámetro de suavizado más bajo, en línea roja se representa en el gráfico la regresión por splines que se obtendría con un parámetro de suavizado de valor \(0.10\).

Aproximación por series de fourier

La forma de Fourier permite aproximar arbitrariamente cerca tanto a la función como a sus derivadas sobre todo el dominio de definición de las mismas. La idea que subyace en este tipo de aproximaciones (que podrían denominarse semi-no-paramétricas) es ampliar el orden de la base de expansión, cuando el tamaño de la muestra aumenta, hasta conseguir la convergencia asintótica de la función aproximante a la verdadera función generadora de los datos y a sus derivadas (Gallant, A.R.; 1981, 1984).

Un polinomio de Fourier viene dado por la expresión:

\(\frac \eta 2 + \sum_{j=1}^k[a_j\cos(jw_0t)+b_j\sin(jw_0t)]\)

donde \(k\) es es el número de ciclos teóricos o armónicos que consideramos, siendo el máximo \(N/2\), \(w_0=\frac {2\pi} N\) es la frecuencia fundamental (también denominada frecuencia angular fundamental). \(t\) toma los valores enteros comprendidos entre \(1\) y \(N\) (es decir, \(t = 1, 2, 3, ...N\)).

Los coeficientes de los armónicos vienen dados por las expresiones:

\[\frac \eta 2 = \frac 2 N \sum_{i=1}^N y_i\], \[a_j=\frac 2 N \sum_{i=1}^N y_i \cos(jw_0t)\],\[b_j=\frac 2 N \sum_{i=1}^N y_i \sin(jw_0t)\]

La aproximación a una función no periódica \(m(x_i;\theta)\) por una serie de expansión de Fourier se se escribe como:

\(m(x_i;\theta)= \eta + \sum_{j=1}^J[a_j\cos(jx_i)+b_j\sin(jx_i)]\)

El vector de parámetros es \(\theta=(\eta, a_1, b_1,...,a_J,b_J)\) de longitud \(K=1+2J\).

Suponiendo que los datos siguieran el modelo \(y_i=m(x_i;\theta)+e_i\), para \(i=1,2,…,N\) estimariamos \(\theta\) por mínimos cuadrados, minimizando:

\[S_n(\theta)=\frac 1 N \sum_{i=1}^N[y_i-m_K(x_i;\theta)]^2\]

Dado que la variable exógena \(x_i\) no esta expresada en forma periódica, debe de transformase o normalizarse en un intervalo de longitud menor que \(2\pi\),\([0,2\pi]\).

Como ejemplo vamos a utilizar la base de datos de la Agencia Española de Meteorológica (Aemet) desde el R-package fda.usc. La base de datos contiene mediciones diarias de temperatura, velocidad del viento y precipitaciones de 73 diferentes estaciones meteorológicas de España para los años 1980 a 2009. En este ejemplo vamos a analizar las temperaturas medias diarias de Santander que representamos gráficamente en R, con la siguiente programación:

library(fda) 

## Loading required package: splines

## Loading required package: Matrix

## 
## Attaching package: 'fda'

## The following object is masked from 'package:graphics':
## 
##     matplot

library(fda.usc)

## Loading required package: mgcv

## Loading required package: nlme

## This is mgcv 1.8-17. For overview type 'help("mgcv-package")'.

data(aemet,package = "fda.usc")
tt = aemet$temp$argvals
temp = as.data.frame(aemet$temp$data,row.names=F)
range.tt = aemet$temp$rangeval
inv.temp = data.frame(t(aemet$temp$data)) # 365 x 73 matrix
names(inv.temp) = aemet$df$name 
plot(ts(inv.temp[,21]),main="Temperaturas medias diarias Santander 1980-2009")

A continuación se van a suavizar estas temperaturas diarias utilizando funciones periódicas de Fourier, en concreto vamos a utilizar las funciones de base igual a 5. Es decir, los armónicos que se obtendrían con:

\(\sum_{j=1}^5[a_j\cos(jw_0t)+b_j\sin(jw_0t)]\)

Santander5 = create.fourier.basis(rangeval = range(tt),nbasis = 5)
plot(Santander5)

La función: smooth.basis(argvals=1:n, y, fdParobj), del R-package fda, donde argvals es el dominio, y es el conjunto de valores a suavizar, y fdParobj, la función base utilizada como regresores:

Santanderfourier5.fd = smooth.basis(argvals = tt, y = inv.temp[,21],fdParobj = Santander5)
plot(ts(inv.temp[,21]),main="Temperaturas medias diarias Santander 1980-2009")
lines(Santanderfourier5.fd,col="red")

Regresion band spectrum

Nerlove (1964) y Granger (1969) fueron los primeros investigadores en aplicar el analisis espectral a las series de tiempo en economía. El uso del analisis espectral requiere un cambio en el modo de ver las series económicas, al pasaser de la perspectiva del tiempo al dominio de la frecuencia. El analisis espectral parte de la supsición de que cuanquier serie \(X_t\), puede ser transformada en ciclos formados con senos u cosenos:

\(X_t=\eta+\sum_{j=1}^N[a_j\cos(2\pi\frac{ft}n)+b_j\sin(2\pi\frac{ft}n)]\) (1)

donde \(\eta\) es la media de la serie, \(a_j\) y \(b_j\) son su amplitud,\(f\) son las frecuencias que del conjunto de las \(n\) observaciones,\(t\) es un indice de tiempo que va de 1 a N, siendo N el numero de periodos para los cuales tenemos observaciones en el conjunto de datos, el cociente \(\frac{ft}n)\) convierte cada valor de \(t\) en escala de tiempo en proporciones de \(2n\) y rango \(j\) desde \(1\) hasta \(n\) siendo \(n=\frac{N}2\) (es decir, 0,5 ciclos por intervalo de tiempo). Las dinámica de las altas frecuencias (los valores más altos de f) corresponden a los ciclos cortos en tanto que la dinámica de la bajas frecuencias (pequeños valores de f) van a corresponder con los ciclos largos. Si nosotros hacemos que \(\frac{ft}n=w\) la ecuación (1) quedaría, asi :

\(X_t=\eta+\sum_{j=1}^N[a_j\cos(\omega_j)+b_j\sin(\omega_j)]\)(2)

El analisis espectral puede utilizarse para identificar y cuantificar en procesos aparentemente aperiodicos, sucesiones de cicos de periodo de corto y largo plazo. Una serie dada \({X_t}\) puede contener diversos ciclos de diferentes frecuencias y amplitudes, y esa combinación de frecuencias y amplitudes de carcter cíclico la hacer aparecer como un serie no periodica e irregular. De hecho la ecuación (2), muestra que cada observación \(t\) de una serie de tiempo, es el resultado sumar los valores en \(t\) que resultan de \(N\) ciclos de diferente longitud y amplitud, a los que habría que añadir si cabe un termino de error.

Una manera practica de pasar desde el dominio del tiempo al dominio de la frecuencia es pre-multiplicando los datos originales por una matriz ortogonal, \(W\), sugerida por Harvey (1978), con el elemento (j,t)th :

\[\begin{equation} w_{jt} = \left\lbrace\begin{array}{ll}\left(\frac{1}T\right) ^\frac{1}2 & \forall j=1\\ \left(\frac{2}T\right) ^\frac{1}2 \cos\left[\frac{\pi j(t-1)}T\right] & \forall j=2,4,6,..\frac{(T-2)}{(T-1)}\\ \left(\frac{2}T\right) ^\frac{1}2 \sin\left[\frac{\pi (j-1)(t-1)}T\right] & \forall j=3,5,7,..\frac{(T-2)}T\\ \left(\frac{1}T\right) ^\frac{1}2 (-1)^{t+1} & \forall j=T\end{array}\right.\end{equation}\]

La matriz \(W\) tiene la ventaja de ser ortogonal por lo que \(WW^T=I\).

Sea \(a\) un vector n x 1 el modelo transformado en el dominio de la frecuencia esta dado por:

\(\hat a= Wa\)

Engle (1974), demostró que una regresión realizada con las series transformadas en el dominio de la frecuencia, Regresión Band Spectrum (RBS), no alteraba los supuestos básicos de la regresión clásica, cuyos estimadores eran Estimadores Lineales Insesgados y Optimos (ELIO).

\[\begin{equation} y=X\beta+u \end{equation}\] (4)

donde \(X\) es una matriz \(n x k\) de observaciones de \(k\) variables independientes, \(\beta\) es un vecto \(k x I\) de parámetros, \(y\) es un vecto \(n x 1\) de observaciones de la variable dependiente, y \(u\) en un vector \(n x I\) de pertubacciones de media cero y varianza constante, \(\sigma^2\).

El modelo se expresaría en el dominio de la frecuencia aplicando una transformación lineal a las variables dependiente e independientes,por ejemplo, premultiplicando todas las variables por la matriz ortogonal \(W\). La técnica de la RBS,consiste en realizar el analisis de regresión en el dominio de la frecuencia pero omitinedo determinadas oscilaciones periodicas. Con este procedimiento pueden tratarse problemas derivados de la estacionalidad de las series o de autocorrelación en los residuos. Engle (1974) muesta que si los residuos están correlacionados serialmente y son generados por un procieso estacionario estocastico, la regresión en el dominio de la frecuencia es el estimador asintóticamente más eficiente de \(\beta\).

La transformación de la ecuación (4) del dominio del tiempo al dominio de la frecuencia en Engle (1974), se basa en la matriz \(W\), cuyo elemento \((t, s)\) esta dado por:

\(w_{ts}=\frac{1}{\sqrt n} e^{i\lambda_t s}\),\(s= 0,1,...,n-1\)

donde \(\lambda_t = 2\pi \frac{t}n\), \(t=0,1,...,n-1\), y \(i=\sqrt{-1}\).

Premultiplicando las observaciones de (4) por \(W\), obtenemos: \[\begin{equation} \dot y=\dot X\beta+\dot u \end{equation}\] (5)

donde \(\dot y = Wy\),\(\dot X = WX\), y \(\dot u = Wu\).

Si el vector de las perturbaciones en (4) cumple las hipoteis clásicas del modelo de regresión: \(E[u] = 0\) y \(E[uu']=\sigma^2 I_n\). entonces el vector de perturbaciones transformado al dominio de la frecuencia, \(\dot u\), tendrá las mimas propiedades. Por otro lado, dado que la matriz W es ortogonal, \(WW^{T}= I\), entonces \(W^T\) sería la transpuesta de la completa conjugada de \(W\). De forma que las observaciones del modelo (5) acaban conteniendo el mismo tipo de información que las observaciones del modelo inicialmente planteado.

\(var(\dot u)=E(\dot u \dot u^T)=E(Wuu'W^T)=WE(uu')W^T=\sigma^2W \Omega W^T\)

si \(\Omega=I\) entonces \(var(\dot u)=\sigma^2I\).

asuminendo que \(x\) es independiente de \(u\), el toerema de Gauss-Markov implicaría que

\(\hat \beta = (\dot x' \dot x)^{-1}\dot x'\dot y\)

es el mejor estimador lineal insesgando (ELIO) de \(\dot \beta\). El estimador obtenido sería de hecho identico al estimador MCO de (4).

Estimar (5) manteniendo unicamente determinadas frecuencias, puede llevarse a cabo omitiendo las observaciones correspondientes a las restantes frecuencias, si bien, dado que las variables en (5) son complejas, Engle (1974) sugiere la transformada inversa de Fourier para recomponer el modelo estimado en términos de tiempo.

Definiendo una matriz de tamaño \(A\) de tamaño n x n de ceros excepto en las posiciones de la diagonal principal correspondientes a las frecuencias que queremos incluir en la regresión y premultiplicando \(\dot y\) y \(\dot X\) por \(A\) eleminamos determindas observaciones y las reemplazamos por ceros para realizar la regresión band spectrum. Devolver al dominio del tiempo estas observaciones requiere:

\[\begin{equation} y^* = W^{T}A\dot y = W^{T}AWy \\ x^* = W^{T}A\dot x = W^{T}AWx \end{equation}\] (6)

Al regresar \(y^*\) sobre \(x^*\) obtenemos un \(\beta\) identico al estimador que obtendríamos al estimar por MCO \(\dot y\) frente a \(\dot x\).

Cuando se realiza la regresión band spectrum de esta manera, ocurre un problema asociado a los grados de libertad de la regresión de \(y^*\) sobre \(x^*\) que asumen los programas estadisticos convencionales, \(n - k\), en vez de los grados de libertad reales que serían unicamente \(n'- k\), donde \(n'\) es el numero de frecuencias incluidas en la regresión band spectrum.

Tan H.B and Ashley R (1999), señalan que el procedimiento de elaboración de una RBS consta de tres etapas:

1.- Transformar los datos originales del dominio del tiempo al dominio de la frecuencia utilizando series finitas de senos y cosenos. Implicaría premultiplicar los datos originales por una matriz ortogonal, \(W\), sugerida por Harvey (1978).

2.- Permitir la variación de a través de m bandas de frecuencia usando variables Dummy \(D_j^1,..D_j^m\). Estas variables se elaboran a partir de submuestras de las \(n\) observaciones del dominio de frecuencias, de esta forma \(D_j^s=\dot x_{j,k}\) si la observación \(j\) está en la banda de frecuencias \(s\) y \(D_j^s=0\), en el resto de los casos.

3.- Re-estimar el resultado del modelo de regresión en el dominio del tiempo con las estimaciones y los coeficientes de las \(m\) variables Dummy. Implicaría premultiplicar la ecuación de regresión ampliada por las variables Dummy por la transpuesta de \(W\).

Realizamos una RBS con datos trimestrales del PIB y empleo en Canada, correspondientes al primer trimestre de 1980 al cuarto trimestre del 2000. En la libreria “descomponer”, las funciones “gdf” y “gdt” transforman las series del tiempo a la frecuencia y de la frecuencia al tiempo, siguiendo la transformación sugerida por Harvey (1978).

library(descomponer)

## Warning: package 'descomponer' was built under R version 3.4.1

## Loading required package: taRifx

library(vars)

## Loading required package: strucchange

## Loading required package: sandwich

## Loading required package: urca

data("Canada")
summary(Canada)

##        e            prod           rw            U        
##  Min.   :929   Min.   :401   Min.   :386   Min.   : 6.70  
##  1st Qu.:935   1st Qu.:405   1st Qu.:424   1st Qu.: 7.78  
##  Median :946   Median :406   Median :444   Median : 9.45  
##  Mean   :944   Mean   :408   Mean   :441   Mean   : 9.32  
##  3rd Qu.:950   3rd Qu.:411   3rd Qu.:461   3rd Qu.:10.61  
##  Max.   :962   Max.   :418   Max.   :470   Max.   :12.77

E <- as.numeric(Canada[, "e"])
PIB <- as.numeric(Canada[, "prod"])-E
summary(lm(E~PIB))

## 
## Call:
## lm(formula = E ~ PIB)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -7.495 -2.363  0.543  2.050  7.123 
## 
## Coefficients:
##             Estimate Std. Error t value             Pr(>|t|)    
## (Intercept) 158.7771    32.6261    4.87            0.0000054 ***
## PIB          -1.4643     0.0608  -24.08 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.24 on 82 degrees of freedom
## Multiple R-squared:  0.876,  Adjusted R-squared:  0.875 
## F-statistic:  580 on 1 and 82 DF,  p-value: <0.0000000000000002

K <- c(rep(1,84))
PIB.1 = gdf(PIB)
E.1=gdf(E)
K.1=gdf(K)
summary(lm(E.1~0+K.1+PIB.1))

## 
## Call:
## lm(formula = E.1 ~ 0 + K.1 + PIB.1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -13.563  -0.784  -0.133   0.420  18.278 
## 
## Coefficients:
##       Estimate Std. Error t value             Pr(>|t|)    
## K.1   158.7771    32.6261    4.87            0.0000054 ***
## PIB.1  -1.4643     0.0608  -24.08 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.24 on 82 degrees of freedom
## Multiple R-squared:     1,   Adjusted R-squared:     1 
## F-statistic: 3.56e+06 on 2 and 82 DF,  p-value: <0.0000000000000002

Teniendo presente el periodograma del PIB y el empleo:

gperiodograma(PIB)

gperiodograma(E)

Utilizamos variables Dummys para separar los ciclos de menor frecuencia del resto, y realizamos una RBS para los ciclos de menor frecuencia.

D1 <- c(rep(1,10),rep(0,74)) 
PIB.1.1=D1*PIB.1
rbs.fit=lm(E.1~0+K.1+PIB.1.1)
plot(ts(E,frequency = 4, start = c(1980, 1)),type="l")
lines(ts(lm(E~PIB)$fitted,frequency = 4, start = c(1980, 1)),type="l",col=2)
lines(ts(gdt(rbs.fit$fitted),frequency = 4, start = c(1980, 1)),type="l",col=3)
legend("top", ncol=3,c("E","Estimado LM","Estimado RBS"),cex=0.6,bty="n",fill=c(1,2,3))

La regresión RBS cuando se utiliza diferenciando o filtrando frecuencias mediante variables Dummys, opera como un filtro frecuencial. Hay que tener presente que las series originales se transformas del tiempo a la frecuencia por la matriz ortogonal que hemos denominado \(W\), y dicha matriz tiene una dimensión (n X n), esto determina que la transformación de los datos de la matriz de regresores para la predicción, ahora de tamaño (n+1 x k), de ser transformados al dominio de la frecuencia por una matriz \(W\) de tamaño (n+1 x n+1), se encontrarían en una frecuencia distinta a la que precisa la separación que ser realiza con las variables Dummys, para solucionar estos problemas de falta de sintonía entre las series de frecuiencia, la predicción ha de realizarse sobre la matriz \(W\) de dimensión (n x n), y las observaciones (2 a n+1) de la matriz de regresores.

Vamos a realizar la predicción correspondiente a un incremento trimestral del empleo del 0,25%. Como la preducción se realiza en el dominio de la frecuencia hay transformar los datos al dominio del tiempo

PIB_2.1980_1.2001=c(PIB[2:84],PIB[84]*1.0025)
D1 <- c(rep(1,10),rep(0,74)) 
PIB.2 = gdf(PIB_2.1980_1.2001)
PIB.2.1=D1*PIB.2
E.2.fitted=rbs.fit$coefficients[1]*K.1+rbs.fit$coefficients[2]*PIB.2.1
E_2.1980_1.2001.fitted=gdt(E.2.fitted)
E_2.1980_1.2001.fitted[84]

## [1] 943.78

plot(ts(E,frequency = 4, start = c(1980, 1)),type="l")
lines(ts(lm(E~PIB)$fitted,frequency = 4, start = c(1980, 1)),type="l",col=2)
lines(ts(E_2.1980_1.2001.fitted,frequency = 4, start = c(1980, 2)),type="l",col=3)
legend("top", ncol=3,c("E","Estimado LM","Estimado RBS"),cex=0.6,bty="n",fill=c(1,2,3))

Regresión en el dominio de la frecuencia con un modelo no lineal.

Asumiendo que las series, \(y_t\),\(x_t\),\(\beta_t\) and \(ut\), pueden ser transformadas en el dominio de la frecuencia:

\[y_t=\eta^y+\sum_{j=1}^N[a^y_j\cos(\omega_j)+b^y_j\sin(\omega_j)\]

\[x_t=\eta^x+\sum_{j=1}^N[a^y_j\cos(\omega_j)+b^y_j\sin(\omega_j)]\]

\[ \beta_t=\eta^\beta+\sum_{j=1}^N[a^\beta_j\cos(\omega_j)+b^\beta_j\sin(\omega_j)]\]

Otenemos dichas series pre-multiplicando (7) por \(W\)

\(\dot y=\dot x\dot\beta+\dot u\) (8)

donde \(\dot y = Wy\),\(\dot x = Wx\), \(\dot \beta = W\beta\) y \(\dot u = Wu\)

El sistema (8) puede reescribirse como:

\[\begin{equation} \dot y=Wx_tI_nW^T\dot \beta + WI_nW^T\dot u \end{equation}\]

Si denominamos \(\dot e=WI_nW^T\dot u\), podrían buscarse los \(\dot \beta\) que minimizaran la suma cuadrática de los errores \(e_t=W^T\dot e\).

Una vez encontrada la solución a dicha optimización, se transformarían las series al dominio del tiempo para obtener el sistema (7).

Para obtener una solución a la minimización de los errores que ofrezca el mismo resultado que la regresión lineal por mínimos cuadrados ordinarios, requiere utilizar una matriz de regresores cuya primera columna sería el vector de tamaño \(N\), \((1,0,0,...)\), la segunda columna sería la primera fila de la matriz \(WI_NX_tW^T\) y las columnas siguientes, corresponderían las a las frecuencias de senos o cósenos que queremos regresar.

Denominando a nueva esta matriz de tamaño (Nxp), ,\(\dot X\), donde \(p=2+j\), siendo \(j\) las frecuencias de seno y coseno elegidas como explicativas, los coeficientes de la solución MCO serían:

\[\dot \beta = (\dot X' \dot X)^{-1}\dot X' \dot y\]

donde \(\beta_{0,1}\) sería el parámetro asociado a la constante,\(\beta_{1,1}\) el asociado a la pendiente, y \(\beta_{1,j}\) los asociados a las frecuencias de senos y cósenos elegidas.

En la libreria descomponer, hay una función “rdf” para realizar dicha regresión.

El algoritmo de calculo “rdf” se realiza en las siguentes fases:

1. Calcula el co-espectro de la serie “x” e “y”

Sea \(x\) un vector n x 1 el modelo transformado en el dominio de la frecuencia esta dado por:

\(\hat x= Wx\)

Sea \(y\) un vector n x 1 el modelo transformado en el dominio de la frecuencia esta dado por:

\(\hat y= Wy\)

Denominando \(p_j\) el ordinal del cross-periodograma de \(\hat x\) y \(\hat y\) en la frecuencia \(\lambda_j=2\pi j/n\), y \(\hat x_j\) el j-th elemento de \(\hat x\) y \(\hat y_j\) el j-th elemento de \(\hat y\), entonces

\[ \left\lbrace \begin{array}{ll} p_j=\hat x_{2j}\hat y_{2j}+\hat x_{2j+1}\hat y_{2j+1} & \forall j = 1,...\frac{n-1}{2}\\ p_j=\hat x_{2j}\hat y_{2j}& \forall j = \frac{n}{2}-1 \end{array} \right . \]

\[p_0=\hat x_{1}\hat y_{1}\]

2. Ordena el co-espectro en base al valor absoluto de \(|p_j|\) y genera un índice en base a ese orden para cada coeficiente de fourier.

3. Calcula la matriz \(Wx_tI_nW^T\) y la ordena en base a el indice anterior.

4. Obtiene \(\dot e=WI_nW^T\dot u\), incluyendo el vector correspondiente al parámetro constante, \((1,0,...0)^n\), y calucula el modelo utilizando los dos primeros regresores de la matriz \(Wx_tI_nW^T\) reordenada y ampliadas, calcula el modelo para los 4 primeros, para los 6 primeros, hasta completar los \(n\) regresores de la matriz.

5. Realiza el Test de Durbin () a los modelos estimados, y selecciona aquellos en donde los \(e_t=W^T\dot e\) están dentro de las bandas elegidas a los niveles de significación que determina la función: \(\alpha=0.1\).

6. De todos ellos elige aquel que tiene menos regresores. Si no encuentra modelo devuelve la solución MCO.

Denominando \(p_j\) el ordinal del periodograma de \(\hat a\) en la frecuencia \(\lambda_j=2\pi j/n\), y \(\hat a_j\) el j-th elemento de \(\hat a\), entonces

\[ \left\lbrace \begin{array}{ll} p_j=\hat a_{2j}^{2}+\hat a_{2j+1}^{2} & \forall j = 1,...\frac{n-1}{2}\\ p_j=\hat a_{2j}^{2}& \forall j = \frac{n}{2}-1 \end{array} \right .\]

\[p_0=\hat a_{1}^{2}\]

Entonces el cuadrado del \(\hat a\) puede ser utilizado como un estimador consistente del periodograma de \(a\).

El test de Durbin esta basado en el siguiente estadistico: \(s_j=\frac{\sum_{r=1} ^j p_r}{\sum_{r=1}^m p_r}\)

donde \(m=\frac{1}{2}n\) para \(n\) par y \(\frac{1}{2}(n-1)\) para \(n\) impar.

El estadístico \(s_j\) ha en encontrarse entre unos límites inferior y superior de valores críticos que han sido tabulados por Durbin (1969). Si bien hay que tener presente que el valor \(p_o\) no se considera en el cálculo del estadístico esto es, \(p_o=\hat v_1=0\)

rdf.mod=rdf(E,PIB)
str(rdf.mod)

## List of 4
##  $ datos      :'data.frame': 84 obs. of  4 variables:
##   ..$ Y  : num [1:84] 930 930 930 931 933 ...
##   ..$ X  : num [1:84] -524 -525 -527 -527 -528 ...
##   ..$ F  : num [1:84] 930 931 932 932 931 ...
##   ..$ res: num [1:84] -0.251 -0.71 -1.344 -0.308 1.368 ...
##  $ Fregresores: num [1:84, 1:18] 1 0 0 0 0 0 0 0 0 0 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:84] "X1" "X2" "X3" "X4" ...
##   .. ..$ : chr [1:18] "C" "1" "2" "3" ...
##  $ Tregresores: num [1:84, 1:18] 0.109 0.109 0.109 0.109 0.109 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:18] "C" "1" "2" "3" ...
##  $ Nregresores: num 18

plot(ts(E,frequency = 4, start = c(1980, 1)),type="l")
lines(ts(lm(E~PIB)$fitted,frequency = 4, start = c(1980, 1)),type="l",col=2)
lines(ts(rdf.mod$datos$F,frequency = 4, start = c(1980, 2)),type="l",col=3)
legend("top", ncol=3,c("E","Estimado LM","Estimado RBS"),cex=0.6,bty="n",fill=c(1,2,3))

Redes Neuronales Artificiales

La RNA es una técnica de inferencia basada en Machine Learnig (aprendizaje máquina). Existen distintos modelos de redes neuronales, los más utilizados el del “Perceptron” y el “backpropagation”.

Hay que tener presente que una red neuronal se considera una “caja negra”, donde lo importante es la predicción, y no cómo se hace.

El proceso incluye como ocurre con este tipo de técnicas una fase de entrenamiento para la optimización de las predicciones, y otra de test o prueba.

Los elementos de la red son:

Las neuronas o nodos
Las capas
    De entrada
    De salida
    Oculta (puede tener a su vez varias capas)
Los pesos
La función de combinación
La función de activación
El objetivo (target)

Los nodos (neuronas) de la capa de entrada, se combinan con los nodos de la capa oculta mediante la función de combinación, que suele ser una combinación lineal de los nodos de entrada mediante los pesos. A las neuronas de las capas ocultas, se aplica una función de activación, que suele ser la tangente hiperbólica de la anterior combinación más un parámetro por nodo oculto, con la que estimamos las neuronas de la capa de salida, y sus errores.

Red Neuronal Artificial

El siguiente script utiliza el package neuralnet para hacer una regresion y predecir el valor de viviendas (expresado en $1000), usando el data set Boston incluido en el package MASS. Para un mejor ajuste del modelo (es decir, para que la red aprenda mejor), se hace un preprocesamiento de los datos, donde se normaliza el datas set Boston. En consecuencia, una vez se obtiene la predicción con la red neural, al presentarse en valores normalizados, hay que transformar los datos para obtener la valoración en términos de precios de la vivienda.

Como es habitual en los ejercicios de minería de datos, el conjunto de datos se divide en una muestra de entrenamiento (el 70%) y otra de test.

La programación procede de : http://apuntes-r.blogspot.com.es/2015/12/regresion-con-red-neuronal.html, en donde se destacan las siguientes notas:

• El parametro threshol = 0.05 indica que las iteraciones se detendran cuando el “Cambio” del error sea menor a 5% entre una iteracion de optimizacion y otra. Este “Cambio” es calculado como la derivada parcial de la funcion de error respecto a los pesos.

• El parametro algorithm = “rprop+” refiere al algoritmo “Resilient Backpropagation”, que actualiza los pesos considerando únicamente el signo del cambio, es decir, si el cambio del error es en aumento (+) o disminución (-) entre una iteración y otra. Para detalles ver: https://en.wikipedia.org/wiki/Rprop

• El parametro hidden = c(7,5) especifica una primera capa oculta con 7 neuronas y una segunda capa oculta con 5 neuronas.

#LIBRERIAS Y DATOS
# -----------------------------------------------------
library(MASS); library(neuralnet); library(ggplot2)

## 
## Attaching package: 'ggplot2'

## The following objects are masked from 'package:psych':
## 
##     %+%, alpha

set.seed(65)
datos    <- Boston
n        <- nrow(datos)
muestra  <- sample(n,n*.70)
train    <- datos[muestra, ]
test     <- datos[-muestra, ]
 
 
# NORMALIZACION DE VARIABLES
# -----------------------------------------------------
maxs      <- apply(train, 2, max)
mins      <- apply(train, 2, min)
datos_nrm <- as.data.frame(scale(datos, center = mins, scale = maxs - mins))
train_nrm <- datos_nrm[muestra, ]
test_nrm  <- datos_nrm[-muestra, ]
 
 
# FORMULA
# -----------------------------------------------------
nms  <- names(train_nrm)
frml <- as.formula(paste("medv ~", paste(nms[!nms %in% "medv"], collapse = " + ")))
 
 
# MODELO
# -----------------------------------------------------
modelo.nn <- neuralnet(frml,
                       data          = train_nrm,
                       hidden        = c(7,5), # ver Notas para detalle 
                       threshold     = 0.05,   # ver Notas para detalle
                       algorithm     = "rprop+" 
                       )
 
 
# PREDICCION
# -----------------------------------------------------
pr.nn   <- compute(modelo.nn,within(test_nrm,rm(medv)))
 
# se transforma el valor escalar al valor nominal original
medv.predict <- pr.nn$net.result*(max(datos$medv)-min(datos$medv))+min(datos$medv)
medv.real    <- (test_nrm$medv)*(max(datos$medv)-min(datos$medv))+min(datos$medv)
 
 
 
# SUMA DE ERROR CUADRATICO
# -----------------------------------------------------
(se.nn <- sum((medv.real - medv.predict)^2)/nrow(test_nrm))

## [1] 8.730475173

#GRAFICOS
# -----------------------------------------------------
# Errores
qplot(x=medv.real, y=medv.predict, geom=c("point","smooth"), method="lm", 
      main=paste("Real Vs Prediccion. Summa de Error Cuadratico=", round(se.nn,2)))

## Warning: Ignoring unknown parameters: method

# Red
plot(modelo.nn)

Red Neuronal Artificial Precio Vivienda

Árbol de regresión

El Arbol de decisión puede ser utilizado para clasificar unidades y tambien para realizar una regresión no paramétricas (Árbol de regresión). Si utilizamos la función “rpart” hay que seleccionar en “method” la elección “ANOVA”. En la opción “control” encontramos parámetros opcionales para controlar el crecimiento de los árboles. Por ejemplo, en control = rpart.control (minsplit = 30, cp = 0.001) se le indica que el número mínimo de observaciones en un nodo sea 30 antes de intentar una división y que una división debe pararse cuando el Coste factor de complejidad sea inferior a un 0.001, con el método anova, esto significa que el R-cuadrado total debe aumentar en 0.001 en cada división.

Utilizando una base de datos “cu.summary” de precios de automoviles, paises en donde han sido fabricados, confiabilidad, consumo de gasolina, y gama del vehiculo, realizamos una estimación con un arbol de regresión.

# Regression Tree Example
library(rpart)

# grow tree
fit <- rpart(Mileage~Price + Country + Reliability + Type,method="anova", data=cu.summary)

printcp(fit) # display the results

## 
## Regression tree:
## rpart(formula = Mileage ~ Price + Country + Reliability + Type, 
##     data = cu.summary, method = "anova")
## 
## Variables actually used in tree construction:
## [1] Price Type 
## 
## Root node error: 1354.5833/60 = 22.576389
## 
## n=60 (57 observations deleted due to missingness)
## 
##            CP nsplit  rel error     xerror        xstd
## 1 0.622885266      0 1.00000000 1.03031850 0.175417788
## 2 0.132060610      1 0.37711473 0.52551930 0.100727441
## 3 0.025440940      2 0.24505412 0.36722405 0.080138910
## 4 0.011603894      3 0.21961318 0.38568677 0.090458117
## 5 0.010000000      4 0.20800929 0.40256940 0.089130973

plotcp(fit) # visualize cross-validation results

summary(fit) # detailed summary of splits

## Call:
## rpart(formula = Mileage ~ Price + Country + Reliability + Type, 
##     data = cu.summary, method = "anova")
##   n=60 (57 observations deleted due to missingness)
## 
##              CP nsplit    rel error       xerror          xstd
## 1 0.62288526607      0 1.0000000000 1.0303184965 0.17541778783
## 2 0.13206061011      1 0.3771147339 0.5255192987 0.10072744084
## 3 0.02544093971      2 0.2450541238 0.3672240504 0.08013890951
## 4 0.01160389411      3 0.2196131841 0.3856867743 0.09045811723
## 5 0.01000000000      4 0.2080092900 0.4025694027 0.08913097265
## 
## Variable importance
##   Price    Type Country 
##      48      42      10 
## 
## Node number 1: 60 observations,    complexity param=0.6228852661
##   mean=24.58333333, MSE=22.57638889 
##   left son=2 (48 obs) right son=3 (12 obs)
##   Primary splits:
##       Price       < 9446.5  to the right, improve=0.6228852661, (0 missing)
##       Type        splits as  LLLRLL,      improve=0.5044405322, (0 missing)
##       Reliability splits as  LLLRR,       improve=0.1263005365, (11 missing)
##       Country     splits as  --LRLRRRLL,  improve=0.1243525220, (0 missing)
##   Surrogate splits:
##       Type    splits as  LLLRLL,     agree=0.950, adj=0.750, (0 split)
##       Country splits as  --LLLLRRLL, agree=0.833, adj=0.167, (0 split)
## 
## Node number 2: 48 observations,    complexity param=0.1320606101
##   mean=22.70833333, MSE=8.498263889 
##   left son=4 (23 obs) right son=5 (25 obs)
##   Primary splits:
##       Type        splits as  RLLRRL,      improve=0.43853834880, (0 missing)
##       Price       < 12154.5 to the right, improve=0.25748498350, (0 missing)
##       Country     splits as  --RRLRL-LL,  improve=0.13345695140, (0 missing)
##       Reliability splits as  LLLRR,       improve=0.01637085564, (10 missing)
##   Surrogate splits:
##       Price   < 12215.5 to the right, agree=0.812, adj=0.609, (0 split)
##       Country splits as  --RRLRL-RL,  agree=0.646, adj=0.261, (0 split)
## 
## Node number 3: 12 observations
##   mean=32.08333333, MSE=8.576388889 
## 
## Node number 4: 23 observations,    complexity param=0.02544093971
##   mean=20.69565217, MSE=2.907372401 
##   left son=8 (10 obs) right son=9 (13 obs)
##   Primary splits:
##       Type    splits as  -LR--L,      improve=0.515359607900, (0 missing)
##       Price   < 14962   to the left,  improve=0.131259377800, (0 missing)
##       Country splits as  ----L-R--R,  improve=0.007022106632, (0 missing)
##   Surrogate splits:
##       Price < 13572   to the right, agree=0.609, adj=0.1, (0 split)
## 
## Node number 5: 25 observations,    complexity param=0.01160389411
##   mean=24.56, MSE=6.4864 
##   left son=10 (14 obs) right son=11 (11 obs)
##   Primary splits:
##       Price       < 11484.5 to the right, improve=0.09693168203, (0 missing)
##       Reliability splits as  LLRRR,       improve=0.07767167054, (4 missing)
##       Type        splits as  L--RR-,      improve=0.04209833909, (0 missing)
##       Country     splits as  --LRRR--LL,  improve=0.02201687475, (0 missing)
##   Surrogate splits:
##       Country splits as  --LLLL--LR, agree=0.80, adj=0.545, (0 split)
##       Type    splits as  L--RL-,     agree=0.64, adj=0.182, (0 split)
## 
## Node number 8: 10 observations
##   mean=19.3, MSE=2.21 
## 
## Node number 9: 13 observations
##   mean=21.76923077, MSE=0.7928994083 
## 
## Node number 10: 14 observations
##   mean=23.85714286, MSE=7.693877551 
## 
## Node number 11: 11 observations
##   mean=25.45454545, MSE=3.520661157

# create additional plots
par(mfrow=c(1,2)) # two plots on one page
rsq.rpart(fit) # visualize cross-validation results  

## 
## Regression tree:
## rpart(formula = Mileage ~ Price + Country + Reliability + Type, 
##     data = cu.summary, method = "anova")
## 
## Variables actually used in tree construction:
## [1] Price Type 
## 
## Root node error: 1354.5833/60 = 22.576389
## 
## n=60 (57 observations deleted due to missingness)
## 
##            CP nsplit  rel error     xerror        xstd
## 1 0.622885266      0 1.00000000 1.03031850 0.175417788
## 2 0.132060610      1 0.37711473 0.52551930 0.100727441
## 3 0.025440940      2 0.24505412 0.36722405 0.080138910
## 4 0.011603894      3 0.21961318 0.38568677 0.090458117
## 5 0.010000000      4 0.20800929 0.40256940 0.089130973

# plot tree
plot(fit, uniform=TRUE,
   main="Regression Tree for Mileage ")
text(fit, use.n=TRUE, all=TRUE, cex=.8)

Bibliografía

Introducción a R: https://www.datacamp.com/courses/introduccion-a-r/?tap_a=5644-dce66f&tap_s=10907-287229

Achim Zeileis, Torsten Hothorn (2002). Diagnostic Checking in Regression Relationships. R News 2 (3), 7-10. URL http://CRAN.R-project.org/doc/Rnews/

Albright,R., Lerman,S. y Manski,C. (1977), “Development Of An Estimation Program For The M. Probit Model”. Federal Highway Administration

Akaike, H. (1974), “A new look at the statistical model identification”, IEEE Transactions on Automatic Control AC-19, pp. 716–723.

Amemiya, T. (1978), “On A Two-Step Estimation Of A Multivariate Logit Model”, Journal Of Econometrics 8.

Ashley, Richard A. (1984), “A Simple Test for Regression Parameter Instability,” Economic Inquiry 22, No. 2, 253-267.

Aznar, A. y Trívez, F. J. (1993), Métodos de Predicción en Economía II: Análisis de Series Temporales, Ed. Ariel.

Bassmann, R. (1957). “A Generalized Classical Method Of Linear Estimation Of Coefficients In A Structural Equation.” Econometrica 25, pp. 77-83

Beltran, Mauricio (2015): “Diseño e implementación de un nuevo clasificador de préstamos bancarios a través de minería de datos”. Tesis Doctoral. Departamento de Economía Aplicada y Estadística. UNED.

Breiman, L. (1996). “Bagging predictors”. Machine Learning 24 (2): 123–140. doi:10.1007/BF00058655. CiteSeerX: 10.1.1.32.9399. http://link.springer.com/article/10.1007%2FBF00058655

Cayuela L (2010) Modelos lineales generalizados (GLM). EcoLab, Centro Andaluz de Medio Ambiente, Universidad de Granada. Junio 2010.

Durbin, J. y Koopman, S. J. (2001), Time Series Analysis by State Space Models (Oxford Statistical Science Series, nº 24), Oxford University Press.

Durbin, J. y Watson, G. S. (1950), “Testing for Serial Correlation Least Squares Regressions”, Biometrika, vol 37. pp. 409-428.

Engle, Robert F. (1974), Band Spectrum Regression,International Economic Review 15,1-11.

Bradley Efron, Elizabeth Halloran, and Susan Holmes (1996). “Bootstrap confidence levels for phylogenetic  trees”. PNAS 93 (23): http://www.pnas.org/content/93/23/13429.full.pdf

Fisher, R. A. (1936). “The Use of Multiple Measurements in Taxonomic Problems”. Annals of Eugenics 7 (2): 179–188.

Fix, E.; J.L. Hodges (1989) “(1951): An Important Contribution to Nonparametric Discriminant Analysis and Density Estimation: Commentary on Fix and Hodges (1951)”. International Statistical Review / Revue Internationale de Statistique 57 (3): 233-238.

Freund, Y; Schapire, R (1997); A Decision-Theoretic Generalization of On-Line Learning and an Application to Boosting, Journal of Computer and System Sciences, 55(1):119-139. http://cseweb.ucsd.edu/~yfreund/papers/adaboost.pdf

Fukunaga y Kessell (1973): “Nonparametric Bayes error estimation using unclassified samples”. IEEE Transactions on Information Theory (Volume:19 , Issue: 4 ):434-440.

Gallant, A. R.(1981) “On the Bias in Flexible Functional Forms and an Essentially Unbiased Form.” J. Econometrics 15(1981):211-45.

Gallant, A. R.(1984) “The Fourier Flexible Form.” Amer. J. Agr. Econ. 66(1984):204-15

Goldfield, S. M. y Quandt, R. E. (1965), “Some test for Homocedasticy”, Journal of American Statistical Association. Vol 37. pp 539-547.

Greene, W. H. (2000), Análisis Econométrico, Ed. Prentice Hall

Gujarati, D. (1997), Basic Econometrics, McGraw-Hill

Gujarati, D. (2003), Econometría, Ed. McGraw-Hill

Hair, J.F., Anderson R.E., Tatham R.L., Black W.C. (2008): Análisis Multivariante. 5ª Edición. Pearson, Prentice Hall.

Hastie, T., Tibshirani R. and Friedman, J. (2008), The Element of Statistical Learning. Data Minining, Inference and Prediction. Second Edition. Springe.

Harvey, A.C. (1978), Linear Regression in the Frequency Domain, International Economic Review, 19, 507-512.

Johnston, J. (1997), Econometric Methods. McGraw-Hill.

Johnston, J. y Dinardo, J. (2001), Métodos De Econometría, Ed. Vicens-Vives 3ª Ed.

Mantegna, R. N.(1997): Degree of Correlation Inside a Financial Market in [Proc. of the ANDM 97 International Conference], Edited by J. Kadtke, AIP press.

Mood, A. M. (1950), Introduction to the Theory of Statistics, McGraw-Hill.

Muñoz A., Parra F. (2007): Econometría Aplicada. Ediciones Académicas

Nelder, John; Wedderburn, Robert (1972). “Generalized Linear Models”. Journal of the Royal Statistical Society. Series A (General) (Blackwell Publishing) 135 (3): 370–384.

Novales, A. (1993), Econometría, 2ª Edición, McGraw-Hill.

Parra, F.(2016): Econometria Aplicada I. https://econometria.files.wordpress.com/2014/11/parra-econometria-aplicada-i1.pdf

Parra, F.(2016): Econometria Aplicada II. https://econometria.files.wordpress.com/2015/01/parra-econometria-aplicada-ii5.pdf.

Pindyck, R. S. y Rubinfield, D. L. (1976), Econometric Models and Economic Forecast, McGraw-Hill.

Pindyck, R. S. y Rubinfield, D. L. (1980), Modelos Econométricos, Ed. Labor.

Pulido, A. (1983), Modelos Econométricos, Ed. Pirámide

Rosenblatt, F. (1958): The Perceptron: A Probabilistic Model for Information Storage and Organization in the Brain, Cornell Aeronautical Laboratory, Psychological Review, v65, No. 6, pp. 386–408. doi:10.1037/h0042519.

Stewart, M. y Wallis, K. (1984), Introducción a la Econometría, Alianza Editorial.

Tan, Hui Boon & Ashley, Richard, 1999. “Detection And Modeling Of Regression Parameter Variation Across Frequencies,” Macroeconomic Dynamics, Cambridge University Press, vol. 3(01), pages 69-83, March.

Venables, W. N. y Ripley, B. D. (2002), Modern Applied Statistics with S. 4ª Ed., Springer.

Werbos, P. (1990): Backpropagation through time: what it does and how to do it. Proceedings of the IEEE, Volume 78, Issue 10, 1550 - 1560, Oct 1990, doi10.1109/5.58337

White, H. (1980), An Heteroskedastic-Consistent Regression with Independent Observation. Econometrica 48, pp. 817-838.