UNIVERSIDAD DE GRANADA FACULTAD DE CIENCIAS DEPARTAMENTO DE ESTADÍSTICA E INVESTIGACIÓN OPERATIVA

MÁSTER EN ESTADÍSTICA APLICADA

ESTIMACIÓN DE UN MODELO ADITIVO NO PARAMÉTRICO

Juan Manuel Parra Murciego

Granada, 7 de Octubre de 2011.

2

TRABAJO DE INVESTIGACIÓN Juan Manuel Parra Murciego

ESTIMACIÓN DE UN MODELO ADITIVO NO PARAMÉTRICO

Memoria de segundo ciclo presentada por D. Juan Manuel Parra Murciego y dirigida por la profesora Dra. María Dolores Martínez Miranda

Vo Bo

María Dolores Martínez Miranda

3

Índice general Prólogo ........................................................................................................................ 5 Capítulo 1. Estimación no paramétrica de modelos de regresión aditivos ............... 7 1.1.Formulación de un modelo de regresión no paramétrico..................................... 7 1.2. El modelo aditivo no paramétrico ....................................................................... 8 1.3. Estimación del modelo aditivo no paramétrico ................................................. 10 Capítulo 2. El método de estimación Backfitting ..................................................... 12 Capítulo 3.El método de estimación Smooth Backfitting ........................................ 20 Capítulo 4. Selección del parámetro de suavizado ................................................... 25 4.1. El parámetro de suavizado en los modelos no paramétricos ............................. 25 4.2. Métodos de selección automática .................................................................... 27 4.3. Selección del parámetro de suavizado Smooth Backfitting ............................... 30 Capítulo 5. Modelización aditiva en la práctica con R. ........................................... 33 5.1 Introducción ...................................................................................................... 33 5.2 Ejemplo con datos reales ................................................................................... 37 5.3 Ejemplo con datos simulados ............................................................................. 53 Bibliografía................................................................................................................ 62

4

Prólogo

Los modelos clásicamente usados en Estadística son paramétricos y para ello se ha de suponer que la muestra de observaciones proviene de una familia paramétrica conocida. En estos casos, el problema es estimar los parámetros desconocidos o hallar tests de hipótesis o intervalos de confianza para los mismos. Esta suposición puede ser relativamente fuerte porque el modelo paramétrico supuesto puede no ser el correcto ya que los datos pueden ser tales que no exista una familia paramétrica adecuada que proporcione un buen ajuste. Por otra parte, los métodos estadísticos desarrollados para un modelo paramétrico particular pueden llevar a conclusiones erróneas cuando se aplican a un modelo ligeramente perturbado (falta de robustez de los datos respecto del modelo). Estos problemas llevaron a la tendencia de desarrollar métodos no paramétricos o semiparamétricos para analizar los datos.

Un modelo paramétrico razonable produce inferencias precisas mientras que un modelo erróneo posiblemente conducirá a conclusiones equivocadas. Por otro lado, los modelos no paramétricos si bien están asociados con alta estabilidad tienen menor precisión. Recientemente, los modelos no paramétricos han ganado una importante atención en el estudio de fenómenos naturales con comportamiento de complejidad no lineal. Sean

el conjunto de pares de datos que provienen de n-observaciones

independientes e idénticamente distribuidas tales que modelo no paramétrico

y

donde los errores

, y sea el son variables

aleatorias independientes e idénticamente distribuidas con media cero, varianza unidad e independientes de las variables aleatorias

. El análisis de este modelo requiere de

técnicas de suavizado multivariadas para la función m y por lo tanto, encuentra a la hora de aplicarlo el problema conocido como “la maldición de la dimensionalidad” que está asociado al hecho de que cuando estamos estimando considerando un entorno con un número fijo de datos, y tenemos una superficie de gran dimensión, dicho entorno puede ser demasiado grande como para ser llamado local; hecho que produce grandes sesgos. Es decir, se necesita un número exponencialmente mayor de datos para que dichos entornos contengan observaciones de la muestra.

5

En los últimos años, para resolver este problema, diversos autores han tratado el problema de reducción de la dimensión de las covariables en modelos de regresión no paramétrica. Hastie y Tibshirani (1990) introdujeron los modelos aditivos que generalizan los modelos lineales, resuelven el problema de “la maldición de la dimensión” y además son de fácil interpretación. Este nuevo planteamiento combina la flexibilidad de los modelos no paramétricos con la simple interpretación del modelo lineal estándar. En el mismo, se supone que

donde las funciones

son las cantidades a estimar. Estimadores para este modelo han sido ampliamente estudiados en la literatura. En nuestro caso nos centraremos en el estudio de los modelos aditivos no paramétricos, especialmente en la técnica de estimación Backfitting y su variante Smooth Bakfitting.

La memoria se estructura en cinco títulos. En el primer título se introduce el concepto de modelo aditivo y su formulación general presentando las ventajas e inconvenientes de este tipo de modelos. Además se describen de forma general los métodos de ajuste más habituales, aunque existen muchos otros. En el segundo capítulo se desarrolla el método de estimación Backfitting, planteándose el problema general que resuelve de forma iterativa así como las fases que componen el algoritmo así como su aplicación práctica mediante polinomios locales. En el tercer capítulo se desarrolla una versión mejorada de método de estimación Backfitting, el Smooth Backfitting, tanto a través del estimador de Nadaraya-Watson como del estimador lineal local. En el cuarto capítulo se formulan diferentes métodos para seleccionar el ancho de banda óptimo para poder aplicar estos estimadores al ajuste no paramétrico de regresión. Por último, en el quinto capítulo, se desarrollan varios ejemplos prácticos tanto con datos reales como simulados de ajuste de modelos aditivos que nos permiten comparar la efectividad de dichos modelos usando diferentes métodos de estimación en el software libre R.

6

Capítulo 1. Estimación no paramétrica de modelos de regresión aditivos

1.1.Formulación de un modelo de regresión no paramétrico

Consideramos un modelo de regresión múltiple no paramétrico habitualmente formulado como:

donde

son observaciones independientes e idénticamente

distribuidas de una variable p+1 dimensional variables aleatorias independientes con

y los residuos

son

para todo i=1…n. Aquí

m es la función de regresión y que modeliza la relación de dependencia entre la variable explicativa X y la variable de respuesta Y.

El objetivo en el problema de regresión es la estimación de la función de regresión m que suponemos desconocida. En dicha tarea la formulación clásica consiste en asumir que m pertenece a alguna clase paramétrica de funciones, con lo cual el problema se resolvería encontrando aquella función dentro de la clase que sea óptima bajo algún criterio de optimalidad (habitualmente el criterio de mínimos cuadrados). En este trabajo estamos interesados en la solución no paramétrica del problema consistente en asumir tan sólo que dicha función es “suave”, entendiendo esta suavidad en términos de derivabilidad.

Para la estimación no paramétrica de la función m, que bajo la formulación inicial se trata de una función p-valuada, existen diversos métodos en la literatura en las últimas décadas. Entre ellos destacamos los métodos basados en núcleos o kernels, como son el estimador de Nadaraya-Watson y en general los métodos basados en estimación polinomial local y los métodos basados en splines. Los primeros se conocen como métodos de suavizado (del inglés smoothing) e involucran en su definición la denominada función núcleo o kernel, denotada por K(t), y habitualmente definida como 7

una densidad simétrica y con soporte compacto; y un parámetro positivo, denotado por h, que define el tamaño de los entorno locales definidos alrededor de cada punto de estimación. Dicho parámetro se denomina ancho de banda, o también parámetro de suavizado, dado que determina la suavidad de la estimación resultante.

Los métodos basados en splines en términos generales calculan la estimación deseada,

, minimizando un criterio de optimalidad penalizado. Habitualmente la

suma de los cuadrados residual más un término que penaliza la falta de suavidad de las funciones

candidatas. Una revisión de tales métodos puede encontrarse por

ejemplo en los libros de Fan y Gijbels (1996) y Hastie y Tibshirani (1990)

Sin embargo los métodos no paramétricos multidimensionales sufren del conocido problema de la dimensionalidad, en inglés “the curse of dimensionality”. En términos generales el problema que plantean estos métodos es que los estimadores locales, que definen entornos de proximidad para la estimación en cada punto, se encuentran que en dimensiones altas los datos en cada entorno son demasiado escasos con lo cual las estimaciones son muy poco precisas.

Este problema que tiene importantes repercusiones desde el punto de vista teórico y práctico hace que los métodos multidimensionales no sean aconsejables en dimensiones elevadas. El problema se soluciona incluyendo alguna restricción o simplificación en la formulación del modelo, dando lugar a modelos de tipo semiparamétrico o bien asumiendo aditividad en las componentes del modelo de regresión.

En este trabajo nos centramos modelos de regresión no paramétricos aditivos. En lo que sigue definiremos dicho modelo y estudiaremos métodos para su estimación.

1.2. El modelo aditivo no paramétrico

Se dice que una función de regresión m satisface un modelo aditivo si tiene la forma

8

con

para todo i=1…n, y además para asegurar la

identificabilidad de las componentes aditivas

suponemos que

para

todo j=1…d.

Obsérvese que asumir un modelo aditivo supone que las covariables tienen efectos separados en la variable respuesta y vienen presentados por las funciones

,

denominadas componentes del modelo aditivo y que suponemos desconocidas.

Bajo la aproximación anterior de nuevo es posible resolver el problema mediante la aproximación paramétrica clásica o bien estimar las funciones

no paramétricamente.

Como ya hemos descrito anteriormente en este trabajo adoptamos la segunda aproximación de modo que tan sólo asumimos cierta suavidad para las funciones componentes del modelo aditivo

.

Obsérvese sin embargo que en cierto modo el modelo aditivo está a medio camino entre el modelo de regresión lineal múltiple paramétrico (que combina aditivamente transformaciones lineales de las variables,

) y el modelo de regresión múltiple no

paramétrico.

Dichos modelos se caracterizan fundamentalmente porque la naturaleza de los efectos de las variables explicativas sobre la variable de respuesta se considera de forma individual. Esto, obviamente, permite ganar en simplicidad y también en interpretabilidad. Es por ello que con este tipo de modelos se pueden resolver los problemas que hasta ahora estaban planteando los modelos multivariantes no paramétricos, principalmente el problema de la falta de dimensionalidad (se producían grandes sesgos en entornos grandes), la falta de interpretabilidad y el excesivo coste computacional (debido al elevado número de operaciones sólo se aplica a suavizadores multidimensionales de dimensión 2 o 3). Ya que al contrario que éstos, los modelos aditivos consideran cada una de las componentes del modelo por separado utilizando suavizadores univariantes y permiten obtener la estimación global como suma de las 9

estimaciones individuales en cada variable, lo que supone poder interpretar fácilmente los efectos de cada variable independiente sobre la dependiente por separado, sin tener en cuenta el resto de variables independientes del modelo, lo que a menudo resulta deseable en contextos como la Economía o la Sanidad.

1.3. Estimación del modelo aditivo no paramétrico Una vez definido el modelo, el paso siguiente consiste en estimarlo a partir de las nobservaciones disponibles. Obsérvese que ). Además, si el parámetro función

(ya que

y todas las funciones

y

fueran conocidas, excepto la

, entonces esta podría estimarse mediante cualquier estimador no

paramétrico univariante (por ejemplo, mediante un ajuste lineal local). Bastaría con aplicar ese estimador al conjunto de datos

, donde

Esta observación llevará a proponer el algoritmo conocido como Backfitting y su posterior mejora, el algoritmo Smooth Backfitting, para estimar este tipo de modelos. Aunque existen muchos otros métodos, todos ellos parten de la idea del algoritmo Backfitting o del método de integración marginal de Linton y Nielsen (1995).

El abanico de técnicas disponibles para estimar no paramétricamente la función de regresión en un modelo aditivo es muy amplio e incluye, entre otras, las siguientes:

1. Algoritmo Backfitting de Hastie y Tibshirani (1990). 2. Método de integración marginal de Linton y Nielsen (1995). 3. Estimador de Hengartner (1996). 4. Estimador eficiente de Linton (1997). 5. Estimador de Kim, Linton y Hengartner (1997, y posteriormente 1999).

10

6. Algoritmo Backfitting con polinomios locales de Opsomer y Ruppert (1997). 7. Smooth Backfitting de Mammen, Linton and Nielsen (1999). 8. Smooth Backfitting de Nielsen y Sperlich (2005).

A continuación pasamos a describir en detalle los que ocuparán nuestra atención en la parte práctica de este trabajo: el algoritmo Backfitting y el algoritmo Smooth Backfitting. En cada uno de ellos se presenta el método general así como sus propiedades más relevantes. Además en la parte práctica de este trabajo se muestran sus aplicaciones tanto con datos reales como simulados con el software R, elegido por ser de licencia libre.

11

Capítulo 2. El método de estimación Backfitting

El algoritmo Backfitting es

un método de Gauss-Seidel para ajustar modelos

aditivos usando residuos parciales de suavización iterativamente, desarrollado por Buja, Hastie and Tibshirani (1989). Partiendo del modelo aditivo antes presentado y teniendo en cuenta que aditivas

representa la matriz de suavizamiento, para estimar las funciones

, se trataría de resolver el sistema de ecuaciones en forma matricial:

A partir de esta idea se desarrolla el algoritmo para estimar las funciones

.

Las fases del algoritmo son: 1. Iniciación:

donde

son estimaciones iniciales cualesquiera (por ejemplo

donde los coeficientes

son los estimados en el modelo de regresión lineal

múltiple.

2. Para cada

donde

obtener:

es un suavizador dado previamente.

12

3. Repetir el paso 2 hasta obtener el valor de convergencia establecido por criterio previamente. El algoritmo involucra la elección de unas funciones iniciales

de las que

en general no se tiene conocimiento a priori y, por ejemplo, se pueden utilizar las funciones de regresión lineales de Y sobre cada una de las variables

. En

cambio, si el algoritmo Backfitting está incluido dentro de otro procedimiento mayor, este es el que proporciona las funciones iniciales. La estimación de forma natural, ya que

, surge

.

A partir de ahora vamos a considerar sin pérdida de generalidad las observaciones, , centradas con respecto a su media, de modo que nos olvidaremos del término constante, α. Nótese que la descripción general del algoritmo Backfitting puede dar la impresión de que es necesario usar el mismo suavizador para cada una de las variables, no obstante, esto se hace así solo para facilitar la presentación, ya que de hecho este método permite que se utilice un suavizador diferente para cada variable, e incluso suavizadores que se apliquen a varias variables o a transformaciones de variables, sin que esto afecte para obtener la convergencia deseada.

Así pues y de forma teórica, siempre que exista la inversa de M es posible escribir la solución de los estimadores directamente mediante la expresión:

Con lo que los suavizadores del modelo aditivo son además suavizadores lineales de las observaciones

. Sea

una matriz por bloques de dimensión

que tiene una

matriz identidad en el bloque d-ésimo y ceros en el resto. La matriz suavizadora aditiva se define como

13

por tanto la estimación para la componente d-ésima del modelo aditivo se puede escribir como

Sea

la matriz de suavizamiento aditiva para la función

completa,

, entonces

Además sea

la matriz de suavizamiento para la función

definida como

variante,

que se puede considerar como el suavizador de un

modelo aditivo en el que los datos se han generado por el siguiente modelo:

Según Hastie y Tibshirani (1990), para alcanzar la convergencia del sistema matricial que determina las funciones aditivas

del modelo aditivo en el caso

bivariante basta con que se cumpla que:

para lo cual las matrices suavizadoras . En el caso de que

deben estar centradas, es decir se produce lo que se conoce como

concurvitiy, que es un efecto análogo a la colinealidad de los modelos lineales y por tanto no deseable.

Para el caso bi-variante las matrices suavizadoras aditivas tendrían según Hastie y Tibshirani (1990) la expresión:

14

siempre que las inversas de las matrices necesarias para su cálculo existan. Para el caso D-variante, según Opsomer (2000), el criterio de convergencia sería:

con un suavizador aditivo para la variable d-ésima:

Como se puede deducir al tratarse de un proceso iterativo el estudio teórico resulta complicado, de hecho sólo existe una definición explícita de los estimadores en el caso de dos variables. Para ello Opsomer y Ruppert (1997) presentan el cálculo de los estimadores empleando polinomios locales. Sean

un conjunto

de variables aleatorias independientes e idénticamente distribuidas con valores en

.

Se supone el siguiente modelo bivariante:

Sean

y

los núcleos equivalentes para la regresión polinómica local en

. En el caso de

donde

y

con j=1,2, estos núcleos equivalentes tienen la siguiente expresión:

y

Para K una función núcleo,

el parámetro ancho de banda y

15

donde pj es el orden del polinomio local que se utiliza para ajustar la componente, representa una matríz de suavizamiento, cuyas filas son los núcleos equivalentes en las observaciones

y

, de forma matricial:

Entonces se define el vector de valores ajustados en las observaciones:

donde

y

y

son las soluciones del siguiente conjunto de ecuaciones de

estimación:

donde

con j=1,2.

Estas ecuaciones se resuelven usando el algoritmo Backfitting, y en este caso particular bivariante, se obtiene la siguiente expresión:

siempre que existan

y

. Para

la expresión del estimador es:

16

Nótese que para garantizar la existencia de estos estimadores se han de cumplir una serie de hipótesis:

1. La función núcleo K es acotada, continua, tiene soporte compacto y la primera derivada tiene un número finito de cambios de signo sobre su soporte. Además, impar y

con j=1,2.

2. Las funciones de densidad

,

y

son continuas, acotadas, tienen soporte

compacto y la primera derivada, de cada una de ellas, tiene un número finito de cambios de signo sobre su soporte. Además,

,

y como condición suficiente pero no necesaria

3.

y

, j=1,2, cuando

Bajo estas tres hipótesis son válidas asintóticamente para todos los elementos matriciales las siguientes aproximaciones:

donde

es una matriz cuyo elemento ij-ésimo es:

Además bajo estas tres hipótesis se cumple que la matriz

es invertible

y se tiene que:

17

cuando

existe, entonces es válida la siguiente aproximación sobre todos

los elementos de las matrices:

Con esto queda demostrado que las matrices inversas anteriormente expuestas estaban bien definidas para polinomios locales de grado estimador

, pero análogamente se demuestra para

y y

cuando

, para el

.

Análogamente, para la extensión del caso de d-variables independientes se parte del modelo:

El modelo se ajustará por regresión polinómica local de grado impar. Sea matriz suavizadora para el polinomio local de grado impar

la

correspondiente a la

variable j-ésima, con j=1…d. Además de la primera hipótesis planteada para el caso bivariante son necesarias las siguientes hipótesis:

1. Las funciones de densidad

con j=1..d son continuas, tienen soporte compacto

y

2.

y

3. Las derivadas

Los estimadores

para cada j.

para todo j, cuando

-ésimas de

.

, con j=1..d existen y son continuas.

se definen como las soluciones a las ecuaciones

normales:

18

y con

. La solución aplicando el algoritmo Backfitting existe y es única, si la

matriz P es invertible. Se considerará esto como otra hipótesis adicional a las anteriores. A partir de todo lo anterior se pueden establecer las expresiones para el cálculo del sesgo y de la varianza condicionados de los estimadores de las componentes aditivas de las variables independientes

con j=1...d e i=1…n siendo:

19

Capítulo 3.El método de estimación Smooth Backfitting

El método conocido como Smooth Backfitting fue desarrollado por Mammen, Linton and Nielsen (1999) e introducido como una mejora del algoritmo Backfitting, demostrado por Nielsen y Sperlich (2005), resulta más eficiente, robusto y sencillo que calcular. Concretamente requiere menor coste computacional en grandes dimensiones y dado que usa una rejilla de puntos para su cálculo de los estimadores resulta mucho más rápido que Backfitting incluso usando estimadores kernel.

Éste parte como una proyección del estimador suavizado de Nadaraya-Watson en las dos dimensiones

. Para obtener el estimador hay que minimizar la suma

suavizada de cuadrados:

dónde K es la función núcleo d-variante e

. La minimización se hace bajo

las restricciones:

siendo

es la estimador marginal del núcleo de la densidad. Así pues el estimador

Nadaraya-Watson resultante queda:

20

donde

es el estimador Nadaraya-Watson marginal normalizado y

estimador del núcleo de la densidad

, siendo

es el

la densidad marginal de

.

Las fases del algoritmo son:

1. Poner r=0 y calcular los pesos suavizados sobre los puntos

2. Obtener:

3. Desde

calcular para todos los puntos

4. Si se alcanza el criterio de convergencia establecido parar, en caso contrario incrementar r en una unidad y repetir el algoritmo desde el paso 3.

El criterio de convergencia habitual es comprobar desde

la desigualdad:

21

con

.

Una alternativa al estimador Smooth Backfitting de Nadaraya-Watson es el estimador Smooth Backfitting lineal local. Éste parte como una proyección del estimador lineal local en el espacio de funciones aditivas, definido como:

donde la minimización se realiza en todas las composiciones de funciones aditivas con

pero además sobre todos los

. Las

fases del algoritmo son:

1. Poner r=0 y calcular los pesos suavizados sobre los puntos

2. Obtener:

22

3. Desde

calcular para todos los puntos

4. Desde

calcular

y

a partir de A, B, C, y D. 23

5. Si se alcanza el criterio de convergencia establecido parar, en caso contrario incrementar r en una unidad y repetir el algoritmo desde el paso 3.

El criterio de convergencia habitual es comprobar que:

24

Capítulo 4. Selección del parámetro de suavizado

4.1. El parámetro de suavizado en los modelos no paramétricos

Los métodos de suavizamiento resultan muy útiles para la inferencia en los modelos de regresión. De hecho, la flexibilidad que caracteriza estos métodos permite descubrir las características que realmente subyacen en los datos disponibles. Sin embargo no siempre es fácil saber qué características observadas en el suavizado están realmente ahí, o cuáles son simplemente un artificio resultante de la variabilidad de los datos.

En este sentido una buena elección del parámetro de suavizado constituye un elemento crucial para obtener buenos resultados en la práctica. La metodología clásica inferencial usando estimadores no paramétricos supone una elección del parámetro de suavizado, lo que habitualmente se hace compensando de forma óptima las componentes de sesgo y de varianza del estimador. Esto, derivada de la minimización de algún criterio de error relacionado con el error cuadrático medio de la estimación resultante. En este caso el parámetro que determina el ancho de banda elegido constituiría

una estimación de lo que se podría definir como parámetro de suavizado

óptimo y por tanto el objetivo es estimar a partir de los datos observados un buen estimador para dicho parámetro óptimo habitualmente desconocido.

Consideremos de nuevo el problema típico de regresión no paramétrica multivariante

Donde m se estimará por un suavizador d-dimensional. En general, este suavizador tiene la expresión:

25

que depende de una función núcleo, K, d-dimensional, y de una matriz H, simétrica de dimensión dxd denominada matriz de ancho de banda. La elección de la función núcleo resulta un problema de poca dificultad resuelto habitualmente utilizando elecciones que infieran buenas propiedades al estudio, mientras que la elección del parámetro de suavizado h, que conforma la matriz H, tiene una importancia crucial en el aspecto y propiedades del estimador de la función de regresión. En la práctica, valores distintos de h pueden producir estimadores completamente distintos.

El parámetro de suavizado controla el equilibrio que el estimador no paramétrico de la función de regresión debe mantener entre el buen ajuste a los datos observados y la capacidad de predecir bien observaciones futuras. Valores pequeños de h dan mucha flexibilidad al estimador y le permiten acercarse a todos los datos observados (cuando h tiende a 0 el estimador acaba por interpolar los datos), pero los errores de predicción asociados serán altos. Hay, por tanto, sobreajuste (overfitting). En el caso de que h tome un tamaño moderado no se ajustaría tan bien a las observaciones (tampoco es necesario, dado que los datos pueden contener ruido aleatorio) pero predeciría mejor. En el otro extremo, si h es demasiado grande, tendremos falta de ajuste (underfitting), como puede ocurrir con los modelos paramétricos globales.

Como se ha comentado anteriormente, buscar el valor adecuado del parámetro de suavizado persigue conseguir el equilibrio entre el sesgo y la varianza del estimador. Para h pequeño el estimador es muy variable (aplicado a muestras distintas provenientes del mismo modelo da resultados muy distintos) y tiene poco sesgo (el promedio de los estimadores obtenidos para muestras distintas es aproximadamente la verdadera función de regresión). Si h es grande ocurre lo contrario.

El parámetro de suavizado puede elegirse de forma manual: comparando los resultados obtenidos para distintos valores de h y eligiendo aquel que, a juicio del investigador, de el resultado más satisfactorio visualmente, o el mas informativo (el que mejor resuma la relación existente entre los datos). Esta forma de proceder está sujeta a la opinión subjetiva del usuario y no puede automatizarse, lo que la hace inviable cuando el número de estimaciones no paramétricas que se han de realizar es grande. Se necesitan, pues, métodos automáticos de selección del parámetro de suavizado. En la

siguiente sección describimos algunos de los métodos más habituales. 26

4.2. Métodos de selección automática Entre los métodos más habituales destacan tres grandes grupos: los métodos tipo plug-in, los métodos basados en validación cruzada y por último los métodos basados en aproximaciones bootstrap. De ellos los más utilizados sin duda son los dos primeros.

1) Métodos tipo Plug-in

La idea que subyace en este tipo de técnicas es que a partir de la expresión optima teórica del parámetro de suavizado, se proponen estimadores de lo desconocido que se incrustarán en dicha expresión teórica. No obstante, aunque estos procedimientos comparten una misma filosofía, podemos hacer una clasificación de los mismos en dos categorías: 

Plug-in asintótico:

Parte de expresiones asintóticas teóricas del error (MSE o MISE) y estima lo desconocido. Y uno de los inconvenientes que tiene es que es necesario realizar hipótesis bastante restrictivas sobre la función de regresión desconocida, lo que limita su campo de aplicación en la práctica. Además, es válido para tamaños muestrales grandes porque en caso de tamaños muestrales pequeños dicho procedimiento no es efectivo. 

Plug-in no asintótico:

Debido al problema que presenta el plug-in asintótico en muestras pequeñas, se toma como punto de partida expresiones desarrolladas de los errores teóricos de naturaleza no asintótica, incrustando en ellas las estimaciones propuestas para los términos desconocidos.

Por otro lado, dependiendo de si utilizan iteraciones para llegar a la solución óptima, se pueden distinguir entre: 27



Plug-in iterativo:

Añaden una sucesión de iteraciones con el fin de conseguir una progresiva mejora del parámetro de suavizado. Su inconveniente es que aumenta el coste computacional. 

Plug-in no iterativo o directo:

Son los que no utilizan iteraciones para la obtención del selector del parámetro.

2) Métodos basados en validación cruzada

Al contrario que en la filosofía plug-in, los métodos de selección automática basados en criterios de validación cruzada no parten de expresiones teóricas del error (exactas o asintóticas), sino que se basan en la minimización de alguna función que aproxime el error del estimador, a partir de los datos. Una medida tal la constituye el Averaged Squared Error, definido como el error cuadrático promediado de las observaciones, dado por:

y una estimación vendría dada por la suma residual de cuadrados:

Ya que esta medida no es insesgada del error ASE se proponen otras funciones de la medida RSS que sí sean insesgadas. Entre ellas destacan el conocido criterio de

28

validación cruzada que propone elegir el parámetro de suavizado h minimizando el criterio definido por:

donde

se define como el estimador del tipo considerado pero calculado

usando todas las observaciones excepto la i-ésima,

. Por ello a este método se

le suele denominar más explícitamente leave-one-out least-squares cross-validation.

Este tipo de procedimientos tiene la ventaja de su facilidad de uso e implementación computacional. Sin embargo también presentan importantes inconvenientes ya que en general dan lugar a estimadores con una alta variabilidad. En las últimas décadas han surgido diversas modificaciones que intentan reducir esta variabilidad sin perder la buenas propiedades del método (véase por ejemplo, el texto de Silverman, 1986, para más detalles sobre este problema)

3) Métodos basados en aproximaciones bootstrap

La metodología bootstrap tiene como propósito ganar información acerca de la distribución de un estimador. Sin embargo en regresión no paramétrica la metodología bootstrap es utilizada fundamentalmente para dos tareas: la primera es la de elegir el parámetro de suavizado o ancho de banda y la segunda es la de construir intervalos de confianza para la curva de regresión.

En el caso de los modelos aditivos tenemos el selector plug-in dado por Opsomer y Ruppert (1998) y Opsomer (2000) aplicados al algoritmo Backfitting de Hastie y Tibshinari (1990). El selector basado en el método de validación cruzada propuesto por Kauerman y Opsomer (¿?)como modificación directa del propuesto inicialmente por Craven y Wahba (1979). Para el algoritmo Smooth Backfitting tenemos un selector basado en penalizaciones de las sumas de cuadrados de los residuos, propuesto por Mammen y Park (2004). Dicho selector será descrito en profundidad en la siguiente sección. Finalmente, cabe destacar el uso del método bootstrap para la selección del 29

parámetro de suavizado; en este caso se pueden encontrar referencias en las publicaciones de Härdle y Bowman (1988) para el estimador de Nadaraya-Watson, Hall (1990) para un estimador sencillo tipo núcleo y Martínez-Miranda (2000) para el estimador lineal local.

4.3. Selección del parámetro de suavizado Smooth Backfitting

El mayor problema de utilizar estimadores Smooth Backfitting multidimensionales radica en la elección del ancho de banda apropiado. Para ello, Mammen y Park (2004) introdujeron selectores de ancho de banda basados en penalizaciones de las sumas de cuadrados de los residuos. En el ambiente de la regresión no paramétrica clásica esto equivale asintóticamente al método de validación cruzada. Si consideramos el estimador Smooth Backfitting de Nadaraya-Watson el algoritmo para seleccionar el ancho de banda óptimo sería:

1. Tomar anchos de banda

arbitrarios, poner r=0 y

y calcular

todos los pesos necesarios para el estimador Smooth Backfitting.

2. Hacer el algoritmo para el estimador Smooth Backfitting para obtener

de

j=1…d, para todos los puntos que se necesiten para calcular:

3. Para j=1…d:

a. Calcular

pero con “0” en la diagonal de la matriz

suavizadora.

30

b. Encontrar

que minimiza la validación cruzada basada en las funciones

k=1…d.

c. Cuando el

sea encontrado, actualizar los pesos correspondientes a

través del algoritmo del estimador Smooth Backfitting. Incrementar j en una unidad.

4. Si el valor de la validación cruzada (cv) ha mejorado, incrementar r en una unidad e ir al paso 2, en caso contrario parar.

En el caso de considerar el estimador Smooth Backfitting linear local el algoritmo para seleccionar el ancho de banda óptimo sería:

1. Tomar anchos de banda

que tiendan a ser suaves y poner k=0 y

y calcular todos los pesos necesarios para el estimador Smooth Backfitting.

2. Hacer el algoritmo para el estimador Smooth Backfitting para obtener para todos los puntos

que se necesiten para calcular:

3. Para j=1…d:

a. Calcular

pero con “0” en la diagonal de la matriz

suavizadora.

b. Encontrar

que minimiza la validación cruzada basada en las funciones

.

c. Cuando el

sea encontrado, actualizar los pesos correspondientes a

través del algoritmo del estimador smooth backfitting. 31

4. Si el valor de la validación cruzada (cv) ha mejorado, incrementar r en una unidad e ir al paso 2, en caso contrario parar.

En la práctica el rendimiento depende de n y d, los suavizadores de las funciones subyacentes, y de la correlación entre las variables explicativas, especialmente en el caso de utilizar el suavizador Nadaraya-Watson para el estimador Smooth Backfitting.

32

Capítulo 5. Modelización aditiva en la práctica con R.

5.1 Introducción

En este capítulo desarrollamos algunas aplicaciones prácticas en R de los modelos expuestos en los capítulos anteriores. Las ilustraciones consideran dos tipos de datos, por un lado datos reales habitualmente utilizados en la literatura de los modelos de regresión no paramétricos aditivos, y por otro lado datos simulados.

El programa R es un entorno de análisis y programación estadístico que forma parte del proyecto de software libre GNU (General Public Licence). R está disponible en la dirección http://www.r-project.org. El proyecto R comenzó en 1995 por un grupo de estadísticos de la universidad de Auckland, dirigidos por Ross Ihaka y Robert Gentleman. R, está basado en el lenguaje de programación S, y está diseñado específicamente para la programación de tareas estadísticas en los años 80 por los Laboratorios Bell AT&T. El lenguaje S se considera un lenguaje de programación estadística de alto nivel orientado a objetos.

Frente a otros lenguajes de programación, R, es sencillo, intuitivo y eficiente ya que se trata de un lenguaje interpretado (a diferencia de otros como Fortran, C++, Visual Basic, etc.). Como programa de análisis estadístico, R-base permite realizar tareas estadísticas sencillas habituales y además permite extensiones que implementan técnicas estadísticas avanzadas. De este modo se cubre las necesidades de cualquier analista, tanto en el ámbito de la estadística profesional como en el de la investigación estadística.

R consta de un sistema base pero la mayoría de las funciones estadísticas vienen agrupadas en distintos packages que se incorporan de forma opcional. Para los métodos de regresión no paramétrica existen funciones disponibles en el package básico stats, no obstante, el uso más adecuado de dichos métodos puede conseguirse a través de funciones incorporadas en varios packages adicionales y actualmente disponibles en la

33

web. Entre estos packages destacan kernSmooth, locpol, np, locfit, loess, sm, lowess, gam y recientemente se ha incluido uno para Smooth Backfitting, el sBF.

Este paquete permite ajustar modelos aditivos de cualquier dimensión mediante el estimador Nadaraya-Watson. También permite calcular el peso de la función núcleo, K. Se permite seleccionar la distancia desde origen así como el tipo de función núcleo. El dominio de la función núcleo es centrado en el origen con lo que a mayor distancia desde el origen menor se hace el peso que se obtiene mediante esta función. Permite introducir la matriz de densidades conjuntas así como la matriz de anchos de banda a utilizar para el cálculo de las estimaciones, y en el caso de no introducirlas, la propia función sBF crea ambas matrices. También permite definir la columna de la matriz de datos que va a ser la variable dependiente del modelo, siendo por defecto la primera columna. Por último se permite especificar el número de iteraciones para el método Backfitting (por defecto 100) y el límite para la convergencia (por defecto el valor 0’0001).

Para calcular la función núcleo, sBF utiliza por defecto el método gaussiano utilizando la matriz identidad, descrito por Gasser, T. and Muller, H.G. (1979), aunque permite utilizar para este cálculo también el método unifrom (g=0), epanechnikov (g=1), biweight (g=2) y triweight (g=3):

“ ”

Para la elección del ancho de banda, h, o ventana, span, el paquete sBF utiliza una función que permite calcular el h óptimo comenzando con una ventana arbitraria (por defecto con el valor 20) y actualizando el valor de h según la fórmula implementada de forma interna en la función. Para el caso unidimensional aplica directamente la expresión a cada columna j de la matriz de datos explicativos (X):

34

y en el caso multidimensional, aplica la función density del paquete stats, que calcula las estimaciones de densidad del núcleo y las almacena en una matriz ancho de banda utilizado

así como el

, desde el valor mínimo al valor máximo de cada una de

las columnas de la matriz de datos explicativos (X) y con el número de puntos equidistantes (m), en la que la densidad se calcula, con un valor de 100 puntos equidistantes por defecto, salvo que se especifique otro valor como argumento de la función sBF.

En el caso de no introducir la matriz de densidades conjuntas (PP), la crea mediante los valores obtenidos de una función como producto de núcleos:

z: rango de filas de 1 a m. s: rango de columnas de 1 al valor devuelto por la función choose(d,2). i,j,v,k: matrices determinadas de manera interna en la función sBF. G: matriz de m-filas y d-columnas de datos generados entre el mínimo y el máximo de la matriz Xj con la función seq del paquete base. m: número de puntos equidistantes, en la que la densidad se calcula (por defecto 100). d: número de columnas de la matriz de datos original menos uno. choose: función que devuelve coeficientes binomiales de sus valores absolutos.

Para un primer cálculo de las funciones univariantes estimadas suavizadas calcula para cada variable (número de columnas de X, d) los m-valores y los almacena en una matriz

utilizando el ancho de banda basado en el cálculo univariante antes

descrito:

Y: matriz con el número de filas del fichero de datos y una columna de la variable dependiente, por defecto tomada de la primera columna del fichero de datos, salvo que se cambie el parámetro que determina la columna dependiente en la función sBF.

35

Para el cálculo de la estimación de la constante

para el modelo estimado se utiliza

la expresión:

Por último, a partir de todo lo anterior y según el número de repeticiones marcado (por defecto 100) realiza el proceso iterativo calculando desde j=1 hasta d:

ó

y comprobando en cada iteración si se obtiene el criterio de convergencia

:

AB es una matriz especial a partir de la combinación de dos matrices, A y B que se generan mediante matrices de densidad conjunta Pp y pP de forma interna en cada iteración.

Y: matriz de densidades univariantes de la matriz P de una fila y m columnas.

36

Existen otros programas disponibles para el análisis de estos modelos. Por ejemplo, en MatLab está disponible todo el software para la estimación según los modelos de Opsomer y Ruppert (1997,1998).

Finalmente, a modo práctico en R, se ilustra el uso de los paquetes locfit y sm para el ajuste de modelos multivariantes y gam y sBF para el ajuste de los modelos no paramétricos aditivos y su comparaciónn mediante su RSS.

5.2 Ejemplo con datos reales Datos del etanol

Consideremos el conjunto de datos obtenidos mediante la experimentación con un motor de un cilindro que utiliza como fuente de alimentación el etanol. Los datos están disponibles en el marco de datos etanol contenido en el package locfit de R. En concreto dentro se recogen 88 observaciones correspondientes a tales repeticiones del experimento con el mismo motor de un cilindro. Las variables observadas son:

Según recomendaciones de Gu (1992) los datos son transformados para el análisis tomando la raíz cuadrada de NOx y realizando una traslación de las covariables al intervalo unidad calculando variables transformadas T.E=(E-0.535)/0.697 y T.C=(C7.5)/10.5.

El objetivo es analizar la influencia que la relación de compresión y la relación de equivalencia tienen en la cantidad de gas que emite el motor. Para ello consideramos un modelo de regresión paramétrico bidimensional donde la variable dependiente es el gas emitido y las covariables son la relación de compresión y la relación de equivalencia.

37

En general y sin asumir aditividad, el modelo considerado se puede describir de la forma:

con

para todo i=1…88.

Bajo tal formulación, la función f no está especificada y deberá ser estimada mediante algún método no paramétrico. En nuestro caso nos centramos en los ajustes de tipo polinomial local introducidos en el capítulo 1, sección 1.1.

Como se ha descrito al inicio de este capítulo, actualmente se dispone en R de varios packages para el ajuste y análisis de tales modelos. No obstante la mayoría de las funciones son sólo aplicables para el caso de una sola covariable. La extensión multivariante actualmente se puede realizar mediante las funciones loess o lowess, contenidas en el package stats, además de la función locfit del package del mismo nombre, locfit. Dichas funciones implementan una versión de la estimación polinomial local propuesta por Cleveland y Devlin (1988) y Loader (1999). Además de tales funciones, el package sm permite, a través de sm.regression, el ajuste de tales modelos multivariantes basados en los métodos desarrrollados en Bowman y Azzalini (1997).

El ancho de banda (h) o en general parámetro de suavizado (span) involucrado en estos estimadores se elegirá en cada caso atendiendo a las facilidades ofrecidas por los correspondientes packages o bien siguiendo recomendaciones de sus autores.

A continuación describimos el uso de tales funciones para el ajuste no paramétrico a los datos del etanol, en este caso un modelo bidimensional. En concreto realizaremos dos ajustes multivariantes, en primer lugar utilizando la función locfit y después la función sm.regresion. En todos los casos describiremos a continuación el código en R que permite el ajuste de los modelos en las observaciones así como su representación gráfica.

1. Ajuste multivariante mediante locfit. Consideramos un ajuste local lineal con un ancho de banda definido fijando el número de observaciones. Para ello utilizamos el 38

argumento “nn=0.5” para definir intervalos de ancho variable pero todos ellos conteniendo al 50% de las observaciones. Consideramos el núcleo de Epanechnikov que puede ser especificado a través del argumento “kern”. El código y los resultados se muestran a continuación:

> library(locfit) Loading required package: akima Loading required package: lattice locfit 1.5-6

2010-01-20

> data(ethanol) > datos > #Transformación > datos[,1] datos[,2] datos[,3] names(datos) locfit.etanol locfit.etanol Call: locfit(formula = T.NOx ~ T.C + T.E, data = datos, alpha = 0.5, deg = 2, kern = "epan")

Number of observations:

88

Family: Gaussian Fitted Degrees of freedom:

10.963

Re > fit.1 RSS.1 RSS.1 [1] 0.01364329 > plot(locfit.etanol, type='persp', theta=30, xlab="T.C", ylab="T.E", + zlab="T.NOx", shade=.75, col="blue", + main="Ajuste del modelo multivariante no paramétrico mediante locfit") 39

> plot(datos[,1], fit.1, xlab="Y observado", ylab="Y predicho") > lines(-2:14,-2:14,col='gray')

Ajuste del modelo multivariante no paramétrico locfit

m2

T.NOx

1.4 1.2 0.6

0.8

1.0

Y predicho

1.6

1.8

m1

0.6

0.8

1.0

1.2

1.4

1.6

1.8

2.0

Y observado

2. Ajuste multivariante utilizando sm.regresion. En este caso consideramos de nuevo ajustes de tipo polinomial local con ancho de banda constante calculado según el criterio de validación cruzada. De nuevo usamos el núcleo Epanechnikov:

> library(sm) Package `sm', version 2.2-4.1 Copyright (C) 1997, 2000, 2005, 2007, 2008, A.W.Bowman & A.Azzalini Type help(sm) for summary information

> matrizbivarianteCE etanol.2 etanol.2$h T.C

T.E

0.898377432 0.003087079 > fit.2 RSS.2 RSS.2 [1] 0.001725721 > sm.regression(matrizbivarianteCE, datos[,1], poly.index=1, method="cv", + xlab="T.C", ylab="T.E", zlab="T.NOx") > plot(datos[,1], fit.2, xlab="Y observado", ylab="Y predicho") >lines(-2:14,-2:14,col='gray')

2.0

T.NOx

1.5 1.0 -0.60 -0.62

E T.

25

-0.64

1.4 1.2 1.0 0.8 0.6

Y predicho

1.6

1.8

2.0

-0.66 10

20 15 T.C

0.6

0.8

1.0

1.2

1.4

1.6

1.8

2.0

Y observado

41

A continuación vamos a asumir el siguiente modelo no paramétrico aditivo:

Las funciones

y

no están especificadas y deberán ser estimadas mediante algún

método no paramétrico aditivo. En nuestro caso vamos a utilizar Backfitting con suavizadores locales (loess) implementado en R mediante el package gam y Smooth Backfitting implementado en R mediante el paquete sBF.

3. Ajuste del modelo aditivo mediante Backfitting con regresión local gam, con un span de valor uno y grado del polinomio dos para ambas covariables:

> library(gam) Loading required package: splines

Attaching package: 'gam'

The following object(s) are masked from 'package:locfit':

gam.slist

> etanol.3 fit.3 RSS.3 RSS.3 [1] 0.01754130 > > plot(datos[,1], fit.3, xlab="Y observado", ylab="Y predicho") > lines(-2:14,-2:14,col='gray') > plot(etanol.3, se=TRUE, ylim = c(-30, 30), residuals=TRUE, ask=TRUE)

42

2.0 1.8 1.6 1.4

Y predicho

1.2 1.0 0.8 0.6 0.6

0.8

1.0

1.2

1.4

1.6

1.8

2.0

30 20 10 0 -10 -30

-20

lo(T.E, span = 1, degree = 2)

20 10 0 -10 -20 -30

lo(T.C, span = 1, degree = 2)

30

Y observado

10

15

20

25

-0.66

-0.65

-0.64

T.C

-0.63

-0.62

-0.61

-0.60

T.E

4. Ajuste del modelo aditivo utilizado Smooth Backfitting utilizando como rejilla los valores generados por la función sBF y con un ancho de ventana de valor dos y rejilla los datos de etanol:

> library(sBF) >etanol.4 > #Estimaciones en las observaciones > m1.hat m2.hat m0.hat > fit.4 RSS.4 RSS.4 [1] 0.2085644 > plot(datos[,1], fit.4, xlab="Y observado", ylab="Y predicho") > plot(datos[,2], m1.hat, xlab="x1", ylab="m1")

2.0 1.0

1.5

Y predicho

2.5

3.0

> plot(datos[,3], m2.hat, xlab="x2", ylab="m2")

0.6

0.8

1.0

1.2

1.4

1.6

1.8

2.0

1.0

-0.4

0.0

-0.2

0.5

m2

0.0

m1

0.2

1.5

2.0

Y observado

10

15

20 x1

25

-0.66

-0.65

-0.64

-0.63

-0.62

-0.61

-0.60

x2

A la vista de los resultados podemos concluir que de los cuatro modelos ajustados el que mejores resultados ofrece por tener un RSS más bajo es el modelo multivariante no paramétrico ajustado mediante sm.regression, aunque locfit y gam ofrecen también RSS aceptables observamos que para sBF dado que es un package muy restrictivo en uso, 44

este conjunto de datos no se ajustan bien, lo que hace que el RSS que genera sea alto comparado con los otros tres modelos. Para ello, se pueden observar los gráficos de datos observados frente a datos predichos por los modelos y ver el ajuste de las predicciones conforme a la realidad. Sabemos de trabajos previos (Cleveland y Devlin, 1988) que estos datos no resultan adecuados para los modelos aditivos ya que existe interacción significativa entre las covariables del modelo.

45

Datos Boston Housing

Consideremos en este caso el conjunto de datos obtenidos de la base de datos Boston Housing. Los datos están disponibles en el marco de datos BostonHousing contenido en el package mlbench de R Sean las variables de estudio:

para 506 vías censadas en la ciudad de Boston. A los datos se le eliminaron 6 outliers, tal y como se recomiendan estudios previos como los de Breiman and Friedman (1985) y Opsomer y Ruppert (1998).

El objetivo es analizar cómo afectan al valor medio de las casas en propiedad ocupadas el número medio de habitaciones, la tasa de impuesto, la proporción entre alumno y profesor por distrito escolar y la proporción de población de baja posición. Para ello consideramos un modelo de regresión no paramétrico de cuatro dimensiones donde la variable dependiente es el valor medio de las casas en propiedad y las covariables son el resto de variables anteriormente descritas. En general y sin asumir aditividad el modelo considerado se describe como:

con

para todo i=1…506.

Bajo tal formulación, la función f no está especificada y deberá ser estimada mediante algún método no paramétrico. En este caso al igual que con el conjunto de datos del etanol nos centramos en los ajustes de tipo polinomial local introducidos en el capítulo 1, sección 1.1.

46

Como se ha descrito al inicio de este capítulo, actualmente se dispone en R de varios packages para el ajuste y análisis de tales modelos. No obstante la mayoría de las funciones son sólo aplicables para el caso de una sola covariable. La extensión multivariante actualmente se puede realizar mediante las funciones loess o lowess, contenidas en el package stats, además de la función locfit del package del mismo nombre, locfit. Dichas funciones implementan una versión de la estimación polinomial local propuesta por Cleveland y Devlin (1988) y Loader (1999). Además de tales funciones, el package sm permite, a través de sm.regression, el ajuste de tales modelos multivariantes basados en los métodos desarrrollados en Bowman y Azzalini (1997).

El ancho de banda (h) o en general parámetro de suavizado (span) involucrado en estos estimadores se elegirá en cada caso atendiendo a las facilidades ofrecidas por los correspondientes packages o bien siguiendo recomendaciones de sus autores.

A continuación describimos el uso de tales funciones para el ajuste no paramétrico a los datos de Boston Housing, en este caso un modelo de cuatro dimensiones. En concreto realizaremos un ajuste multivariante utilizando la función locfit ya que en este caso no podemos usar la función sm.regression por estar limitado su uso a un máximo de dos covariables. En todos los casos describiremos a continuación el código en R que permite el ajuste de los modelos en las observaciones así como su representación gráfica.

1. Ajuste multivariante mediante locfit. Consideramos un ajuste local lineal con un ancho de banda definido fijando el número de observaciones. Para ello utilizamos el argumento “nn=0.5” para definir intervalos de ancho variable pero todos ellos conteniendo al 50% de las observaciones. Consideramos el núcleo de Epanechnikov que puede ser especificado a través del argumento “kern”.

>#Cargamos el conjunto de datos sin outliers >boston library(locfit) Loading required package: akima Loading required package: lattice 47

locfit 1.5-6

2010-01-20

> locfit.boston locfit.boston Call: locfit(formula = MV ~ rm + tax + ptratio + lstat, data = boston, alpha = 0.5, deg = 2, kern = "epan")

Number of observations:

500

Family: Gaussian Fitted Degrees of freedom: Residual scale:

48.074

3.14

> fit.1 RSS.1 RSS.1 [1] 7.468223 > > plot(boston[,1], fit.1, xlab="Y observado", ylab="Y predicho")

30 10

20

Y predicho

40

50

> lines(-2:114,-2:114,col='gray')

10

20

30

40

50

Y observado

A continuación vamos a asumir el siguiente modelo no paramétrico aditivo:

48

Las funciones

,

,

y

no están especificadas y deberán ser estimadas

mediante algún método no paramétrico aditivo. En nuestro caso vamos a utilizar Backfitting con suavizadores locales (loess) implementado en R mediante el package gam y Smooth Backfitting implementado en R mediante el paquete sBF.

2. Ajuste del modelo aditivo mediante Backfitting con regresión local gam, con un span de valor uno y grado del polinomio dos para ambas covariables:

> library(gam) Loading required package: splines

Attaching package: 'gam'

The following object(s) are masked from 'package:locfit':

gam.slist

> boston.2 fit.2 RSS.2 RSS.2 [1] 10.44345 > > plot(boston[,1], fit.2, xlab="Y observado", ylab="Y predicho") >lines(-2:114,-2:114,col='gray') > plot(boston.2, se=TRUE, ylim = c(-30, 30), residuals=TRUE, ask=TRUE)

49

50 40 10

20

30

Y predicho

10

20

30

40

50

30 20 10 0 -10 -30

-20

lo(tax, span = 1, degree = 2)

20 10 0 -10 -20 -30

lo(rm, span = 1, degree = 2)

30

Y observado

4

5

6

7

8

200

300

500

600

700

30 20 10 0 -10 -30

-20

lo(lstat, span = 1, degree = 2)

20 10 0 -10 -20 -30

lo(ptratio, span = 1, degree = 2)

400 tax

30

rm

14

16

18 ptratio

20

22

10

20

30

lstat

3. Ajuste del modelo aditivo utilizado Smooth Backfitting utilizando como rejilla los 50 primeros valores del conjunto de datos BostonHousting sin outliers:

50

> library(sBF) >boston.3 > #Estimaciones en las observaciones > m1.hat m2.hat m3.hat m4.hat m0.hat > fit.3 RSS.3 RSS.3 [1] 29152.81 > plot(bostonreducido[,2], m1.hat, xlab="x1", ylab="m1") > plot(bostonreducido[,3], m2.hat, xlab="x2", ylab="m2") > plot(bostonreducido[,4], m3.hat, xlab="x3", ylab="m3")

0

Y predicho

100

200

300

> plot(bostonreducido[,5], m4.hat, xlab="x4", ylab="m4")

15

20

25

30

35

Y observado

51

400

150

m2

200

100 -50

-200

0

0

50

m1

5.5

6.0

6.5

7.0

220

240

260

150 m4

50

100

200 100

-50

-200

0

-100

0

m3

300

x2

300

x1

280

15

16

17

18 x3

19

20

21

5

10

15

20

25

30

x4

A la vista de los resultados podemos concluir que de los cuatro modelos ajustados el que mejores resultados ofrece por tener un RSS más bajo es el modelo multivariante no paramétrico ajustado con locfit. De nuevo observando los gráficos y los RSS que son altos en todos los modelos nos hace pensar la existencia de interacciones entre las covariables. Además los modelos ajustados mediante aditividad de nuevo no funcionan adecuadamente y es posible que sea también por el elevado número de observaciones a la vista de los gráficos. En el caso de sBF el ajuste individual de cada covariable con respecto a las observaciones es muy malo, el modelo presenta problemas de multicolinealidad.

52

5.3 Ejemplo con datos simulados En este caso se va realizar un estudio comparativo de los estimadores, el polinomial local multivariante y el aditivo mediante Backfitting y Smooth Backfitting. Se considera de tamaño muestral n=500 para el cual se generan datos a partir de las funciones propuestas por Opsomer y Ruppert (1998) para la exploración del comportamiento de selector plug-in:

Con lo que le modelo no paramétrico bivariante no paramétrico queda definido de la forma:

donde los residuos covariables (

se han generado mediante una distribución según una

y el vector de

.

> #Simulamos los datos para las variables X_1 y X_2 con n=500 de una > #distribución N(0,1) y lo mismo para el error > > x.grid x.grid x.grid datos_X1 datos_X2 e > Y = (datos_X1)^2 + sin(datos_X2) + e

53

Bajo tal formulación, la función m no está especificada y deberá ser estimada mediante algún método no paramétrico. En este caso al igual que con el conjunto de datos del etanol nos centramos en los ajustes de tipo polinomial local introducidos en el capítulo 1, sección 1.1.

Como se ha descrito al inicio de este capítulo, actualmente se dispone en R de varios packages para el ajuste y análisis de tales modelos. No obstante la mayoría de las funciones son sólo aplicables para el caso de una sola covariable. La extensión multivariante actualmente se puede realizar mediante las funciones loess o lowess, contenidas en el package stats, además de la función locfit del package del mismo nombre, locfit. Dichas funciones implementan una versión de la estimación polinomial local propuesta por Cleveland y Devlin (1988) y Loader (1999). Además de tales funciones, el package sm permite, a través de sm.regression, el ajuste de tales modelos multivariantes basados en los métodos desarrrollados en Bowman y Azzalini (1997).

El ancho de banda (h) o en general parámetro de suavizado (span) involucrado en estos estimadores se elegirá en cada caso atendiendo a las facilidades ofrecidas por los correspondientes packages o bien siguiendo recomendaciones de sus autores.

A continuación describimos el uso de tales funciones para el ajuste no paramétrico a los datos que hemos simulado, en este caso un modelo bidimensional. En concreto realizaremos dos ajustes multivariantes, en primer lugar utilizando la función locfit y después la función sm.regresion. En todos los casos describiremos a continuación el código en R que permite el ajuste de los modelos en las observaciones así como su representación gráfica.

1. Ajuste multivariante mediante locfit. Consideramos un ajuste local lineal con un ancho de banda definido fijando el número de observaciones. Para ello utilizamos el argumento “nn=0.7” para definir intervalos de ancho variable pero todos ellos conteniendo al 70% de las observaciones. Consideramos el núcleo de Epanechnikov que puede ser especificado a través del argumento “kern”.

> library(locfit) Loading required package: akima 54

Loading required package: lattice locfit 1.5-6

2010-01-20

> locfit.simulados locfit.simulados Call: locfit(formula = Y ~ datos_X1 + datos_X2, alpha = 0.5, deg = 1, kern = "epan")

Number of observations:

100

Family: Gaussian Fitted Degrees of freedom: Residual scale:

5.71

0.09

> fit.1 RSS.1 RSS.1 [1] 0.007404376 > plot(locfit.simulados, type='persp', theta=30, xlab="X1", ylab="X2", + zlab="Y", shade=.75, col="blue", + main="Ajuste del modelo multivariante no paramétrico mediante locfit") > plot(Y, fit.1, xlab="Y observado", ylab="Y predicho") > lines(-2:14,-2:14,col='gray')

Ajuste del modelo multivariante no paramétrico locfit

m2

Y m1

55

1.5 1.0

Y predicho

0.5 0.0 0.0

0.5

1.0

1.5

Y observado

2. Ajuste multivariante utilizando sm.regresion. En este caso consideramos de nuevo ajustes de tipo polinomial local con ancho de banda constante calculado según el criterio de validación cruzada. De nuevo usamos el núcleo Epanechnikov:

> library(sm) Package `sm', version 2.2-4.1 Copyright (C) 1997, 2000, 2005, 2007, 2008, A.W.Bowman & A.Azzalini Type help(sm) for summary information

> matrizbivarianteX1X2 sm.simulados sm.simulados$h datos_X1 datos_X2 0.1722281 0.1722281 > fit.2 RSS.2 RSS.2 [1] 0.006902935 > sm.regression(matrizbivarianteX1X2, Y, poly.index=1, method="cv", + xlab="X1", ylab="X2", zlab="Y") > plot(Y, fit.2, xlab="Y observado", ylab="Y predicho") > lines(-2:14,-2:14,col='gray') 56

1.5

Y

1.0 0.5

0.8 0.6

X2

0.8

0.4 0.2

1.0 0.0

0.5

Y predicho

1.5

0.2

0.6 0.4 X1

0.0

0.5

1.0

1.5

Y observado

A continuación vamos a asumir el siguiente modelo no paramétrico aditivo:

Las funciones

y

no están especificadas y deberán ser estimadas mediante algún

método no paramétrico aditivo. En nuestro caso vamos a utilizar Backfitting con suavizadores locales (loess) implementado en R mediante el package gam y Smooth Backfitting implementado en R mediante el paquete sBF.

3. Ajuste del modelo aditivo mediante Backfitting con regresión local gam, con un span de valor uno y grado del polinomio dos para ambas covariables:

57

> library(gam) Loading required package: splines

Attaching package: 'gam'

The following object(s) are masked from 'package:locfit':

gam.slist

> datos names(datos) dat simulados.gam fit.3 RSS.3 RSS.3 [1] 0.006751879 > > plot(datos[,1], fit.3, xlab="Y observado", ylab="Y predicho") > lines(-2:14,-2:14,col='gray')

1.0 0.5

Y predicho

1.5

>plot(simulados.gam, se=TRUE, ylim = c(-30, 30), residuals=TRUE, ask=TRUE)

0.0

0.5

1.0

1.5

Y observado

58

30 20 10 0 -10 -30

-20

lo(datos_X2, span = 0.5, degree = 1)

30 20 10 0 -10 -20 -30

lo(datos_X1, span = 0.5, degree = 1)

0.2

0.4

0.6

0.8

0.2

0.4

datos_X1

0.6

0.8

datos_X2

4. Ajuste del modelo aditivo utilizado Smooth Backfitting utilizando como rejilla los valores generados por la función sBF y con un ancho de ventana de valor diez, con rejilla los datos simulados:

> library(sBF) > simulados.sBF #Estimaciones en las observaciones > m1.hat m2.hat m0.hat > fit.4 RSS.4 RSS.4 [1] 0.02860518 > plot(datos[,1], fit.4, xlab="Y observado", ylab="Y predicho") > lines(-2:14,-2:14,col='gray') > plot(datos[,2], m1.hat, xlab="x1", ylab="m1") >lines(0:1,-0.5:1,col='gray') > plot(datos[,3], m2.hat, xlab="x2", ylab="m2") > lines(0:1,-0.5:1,col='gray') 59

1.5 1.0 0.5

Y predicho

0.0

0.5

1.0

1.5

0.0

m2

0.2 -0.4

-0.4

-0.2

-0.2

0.0

m1

0.4

0.2

0.6

0.4

Y observado

0.2

0.4

0.6 x1

0.8

0.2

0.4

0.6

0.8

x2

A la vista de los resultados podemos concluir que de los cuatro modelos ajustados el que mejores resultados ofrece por tener un RSS más bajo es el modelo multivariante no paramétrico ajustado con sm.regression. De nuevo observando los gráficos y los RSS que son bajos en todos los modelos nos hace pensar la no existencia de interacciones entre las covariables. Además los modelos ajustados mediante aditividad se ajustan bien para este tipo de datos. Dado que Smooth Backfitting usa aproximaciones constantes (Nadaraya-Watson) mientras que con las otras funciones que usan estimaciones locales de forma estricta no sería posible comparar los resultados que ofrece con el resto de modelos, pero se adjunta a modo ilustrativo como una alternativa al resto de estimaciones locales.

A continuación se muestra una tabla comparativa de los 4 modelos utilizados con sus características principales: 60

Modelo

Estimador

Multivariante

Polinomios

locfit

locales

Multivariante

Polinomios

sm.regression

locales

Aditivo

Polinomios

Backfitting

locales

Aditivo Smooth

Nadaraya-

Backfitting

Watson

K

h

RSS

Vecino más cercano Epanechnikov

para un % fijo de los

0.007404376

datos Epanechnikov

Criterio de Validación Cruzada

No usa núcleo,

Usa tamaño de

utiliza el grado

ventana prefijado, no

del poninomio

ancho de banda

Gaussiano

Penalizaciones en las RSS calculadas

0.006902935

0.006751879

0.02860518

K: función núcleo. h: Ancho de banda. RSS: Suma de los cuadrados de los residuos.

61

Bibliografía [1]

Buja, A., Hastie, T.J. and Tibshirani, R. (1989). Linear smoothers and additive models (with discussion). The Annals of Statistics, Vol. 17, No. 2, pp.453–555.

[2]

Clark, R.M. (1977) Non-parametric estimation of a smooth regression function. Journal of the Royal Statistical Society, Ser. B, Vol. 39, pp. 107-113.

[3]

Cleveland, W.S. (1979). Robust Locally Weighted Regression and Smoothing Scatterplots. Journal of the American Statistical Association, Vol. 74, No. 368. Theory and Methods, pp. 829-836.

[4]

Cleveland, W.S. and Devlin, S.J. (1988) Locally weigthed regression: An approach to regression analysis by local fitting. Journal of the American Statistical Association, Vol. 83, No. 403, pp. 596-610.

[5]

Cleveland, W. and Grosse, E. (1991). Computational Methods for Local Regression. Statistics and Computing 1.

[6]

Fan, J. (1992) Design-adaptive nonparametric regression. Journal of the American Statistical Association, No. 87, pp. 998-1004.

[7]

Fan, J. and Gijbels, I. (1996) Local polynomial modelling and its applications. Chapman & Hall.

[8]

Fan, J., Gijbels, I. and Hu, T.-C. and Huang, L.-S. (1996). An asymptotic study of variable bandwidth selection for local polynomial regression with application to density estimation. Statistica Sinica, Vol. 6, No. 1.

[9]

Fan, J., Härdle, W. and Mammen, E. (1998) Direct estimation of lowdimensional components in additive models. The Annals of Statistics, Vol. 26,No. 3, pp. 943971.

[10] Gasser, T. and Müller, H.G. (1979) Kernel estimation of regression functions.In Smoothing Techniques for Curve Estimation, Lecture Notes in Mathematics, No. 757, pp. 23-68. Springer-Verlag, New York. [11] Härdle W. (1990). Smoothing techniques. Springer Series in Statistics, New York (1991). 62

[12] Hastie, T.J. and Tibshirani, R. (1986) Generalized additive models (with discussion). Journal of the Royal Statistical Society, No. 1, pp. 297-318. [13] Hastie, T.J. and Tibshirani, R. (1987) Generalized additive models: Some applications. Journal of the American Statistical Association, Vol. 82, No. 398,pp. 371-386. [14] Hastie, T.J. and Tibshirani, R. (1990) Generalized additive models. Washington, D.C.; Chapman Hall. [15] John Fox. (2002) Nonparametric Regression. Appendix to An R and S-PLUS Companion to Applied Regression. [16] Jones, M.C., Davies, S.J. and Park, B.U. (1994) Versions of kernel-type regression estimators. Journal of the American Statistical Association, Vol. 89, No. 427, pp. 825-832. [17] Kauermann, G., and Opsomer, J.D. (2003). Local Likelihood Estimation in Generalized Additive Models. Scandinavian Journal of Statistics, 30, 317–337. [18] Linton, O.B. and Härdle, W. (1996) Estimating additive regression models with known links. Biometrika, No. 83, pp. 529-540. [19] Linton, O.B. (1997) Efficient estimation of additive nonparametric regression models. Biometrika, Vol. 84, No. 2, pp. 469-473. [20] Loader, C. (1999). Local Regression and Likelihood. Springer, New York. [21] Mammen, E., Linton, O.B. and Nielsen, J.P. (1999) The existence and asymptotic properties of a backfitting projection algorithm under weak conditions. The Annals of Statistics, Vol. 27, No. 5, pp. 1443-1490. [22] Mammen, E., and Park, C. (2005). Bandwidth Selection for Smooth Back-fitting in Additive Models. The Annals of Statistics, 33, 1260–1294. [23] Martínez Miranda, M.D., Raya Miranda, R., González Manteiga, W. and González Carmona, A. (2008) “A bootstrap local bandwidth selector for additive models”, Journal of Computational and Graphical Statistics, 17, 38-55.

63

[24] Nadaraya, E.A. (1964) On estimating regression. Theory Probab. Appl, No. 9, pp. 141-142. [25] Nielsen, J.P., and Sperlich, S. (2005). Smooth Backfitting in Practise. Journal of the Royal Statistical Society, Ser. B, 67, 43–61. [26] Opsomer, J.D. and Ruppert, D. (1997) Fitting a bivariate additive model by local polynomial regression. The Annals of Statistics, Vol. 25, No. 1, pp.186-211. [27] Opsomer, J.D., and Ruppert, D. (1997). Fitting a Bivariate Additive Model by Local Polynomial Regression. The Annals of Statistics, 25, 186–211. [28] Opsomer, J.D. and Ruppert, D. (1998) A fully automated bandwidth selection method for fitting additive models. Journal of the American Statistical Association, Vol. 93, No. 442, pp. 605-619. [29] Opsomer, J.D. (2000) Asymptotic properties of backfitting estimators. Journal of Multivariate Analysis, Vol. 73, pp. 166-179. [30] Ruppert, D. and Wand, M.P. (1994). Multivariate Locally Weighted Least Squares Regression. The Annals of Statistics, Vol. 22, No. 3, pp. 1346-1370. [31] Ruppert, D., Sheather, S. J. and Wand, M. P. (1995). An effective band- width selector for local least squares regression. Journal of the American Statistical Association, 90, 1257–1270. [32] Ruppert, D. Wand, M.P. and Carroll, R.J. (2003). Semiparametric Regres-sion. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press. [33] Sheather, S. J. and Jones, M. C. (1991). A reliable data-based bandwidth selection method for kernel density estimation. Journal of the Royal Statistical Society, Series B, 53, 683–690. [34] Sheather, S. J. and Jones, M. C. (1991). A reliable data-based bandwidth selection method for kernel density estimation. JRSS-B 53, 683-690. [35] Silverman, B.W. (1986) Density Estimation for Statistics and Data Analysis. Monographs on Statistics and Applied Probability. Ed. Chapman and Hall. 64

[36] Wand, M. P. (1994). Fast Computation of Multivariate Kernel Estimators. Journal of Computational and Graphical Statistics, 3, 433-445. [37] Wand, M.P. and Jones, M.C. (1994). Multivariate Plug-in Bandwidth Selection. Computational Statistics, 9. pp. 97-116. [38] Wand, M.P. and Jones, M.C. (1995). Kernel Smoothing. Monographs on Statistics and Applied Probability. Ed. Chapman and Hall. [39] Yang, L. and Tschernig, R. (1999) Multivariate bandwidth selection for local linear regression. Journal of the Royal Statistical Society, Ser. B, Stat Methodol. 61, No. 4, pp. 793-815.

65