Los siguiente temas especiales han sido desarrollado por estudiantes como parte de su primer proyecto de curso.
Equipo 1:
Iván García
Samuel Gurrola
Rubí Lara
Equipo 2:
Rodolfo Ramírez
Alfredo Bistrain
Jesús Salazar
Origen estadístico del cuantil
Regresión cuantíl: nacimiento formal (1978)
Propuesta formal por Roger Koenker y Gilbert Bassett:
Artículo fundacional: Regression Quantiles, Econometrica, 1978.
Generalización de la regresión de la mediana a cualquier cuantil
\[
\tau \in (0,1)
\]
Utiliza una función de pérdida asimétrica y se resuelve con programación lineal.
Difusión y aplicaciones
Difusión en estadística y econometría
quantreg desarrollado por Koenker en R.Limitación del modelo lineal tradicional
Nota: Bajo el modelo clásico
\[ Y_i = X_i^\top \beta + \varepsilon_i, \quad \mathbb{E}[\varepsilon_i \mid X_i] = 0 \]
se tiene
\[ \mathbb{E}[Y_i \mid X_i] = X_i^\top \beta. \]
Ventaja de la regresión cuantíl
Comparación visual
Fundamentos de la regresión cuantíl
Definición de cuantiles
Conceptos Básicos de la Regresión Cuantílica El cuantil como problema de minimización
En este apartado veremos cómo los cuantiles pueden ser expresados como un problema de minimización. Es bien sabido que la media de \(Y\) es la solución al problema de minimización siguiente: \[ \arg\min_c E[(Y - c)^2] \] Demostraremos más adelante que si reemplazamos la función de pérdida cuadrática por la función de pérdida absoluta, se tiene que la mediana, que denotaremos por \(Q_{Y}\left(\frac{1}{2}\right)\) es solución de:
\[ \arg\min_c E[|Y - c|] \]
Para generalizar el hecho anterior para cualquier cuantil, definimosla función de pérdida: \[ \rho_{\theta}(y) = [\theta - I(y < 0)] \cdot |y| = y(1 - \theta) I(y \leq 0) + y \theta I(y > 0) \] Así, como demostraremos más adelnate, el cuantíl \(Q_{Y}(\theta)\) es solución de la siguiente función de minimización
\[ \arg\min_c E[\rho_{\theta}(Y - c)] \] Ahora demostraremos que la mediana minimiza la esperanza de la función de pérdida del valor absoluto, esto para el caso en el que \(Y\) es una v.a. con densidad \(f(y)\). Así, nótese que,
\[ \begin{align} \mathbb{E}[|Y-c|] &= \int_{y \in \mathbb{R}}|y-c|f(y)dy\\ &= \int_{y<c}(c-y)f(y)dy + \int_{y>c}(y-c)f(y)dy \end{align} \] Analizando la integral de la izquierda de la última expresión se tiene que:
\[ \int_{y<c}(c-y)f(y)dy = c \int_{-\infty}^{c} f(y)dy - \int_{-\infty}^{c} y f(y)dy \] Dado que el valor absoluto es una función convexa entonces derivando respecto a \(c\) e igualando a 0, dbeemos ser capaces de encontrar el punto en el que se encuentra el mínimo. Así,
\[ \begin{align*} \frac{\partial}{\partial c}\left[\int_{y<c}(c-y)f(y)dy\right] &= cf(c)+F(c)-cf(c)\\ &= F(c) \end{align*} \]
De forma similar:
\[ \begin{equation*} \frac{\partial}{\partial c}\left[\int_{y<c}(c-y)f(y)dy\right] = F(c) -1 \end{equation*} \]
Así, juntando los dos resultados anteriores, se tiene que
\[ \frac{\partial}{\partial c}\mathbb{E}[|Y-c|] = 2F(c)-1 \]
por lo que
\[ \frac{\partial }{\partial c}\mathbb{E}[|Y-c|] = 0 \Longleftrightarrow F(c) = \frac{1}{2} \Rightarrow c = Q_{Y}(0.5) \]
Para el caso general, dónde la función de pérdida es \(\rho_{\theta}\) se tiene que :
\[ \mathbb{E}[\rho_{\theta}(Y-c)] = (1- \theta)\int_{y<c}(c-y)f(y)dy + \theta\int_{y>c}(y-c)f(y)dy \]
y podemos ver que es prácticamente la misma expresión a la que llegamos en el caso de la mediana excepto que estas integrales están multiplicadas por unas constantes. Así, bajo el mismo procedimiento enterior, se tiene que \[ \frac{\partial }{\partial c}\mathbb{E}[\rho_{\theta}(Y-c)] = (1- \theta)F(c) + \theta(F(c)-1) = F(c) - \theta \]
por lo tanto, dado que \(\rho_{\theta}\) es una función convexa entonces alcanza su mínimo en dódne la primera derivada se hace 0 y por lo tanto,
\[ \frac{\partial }{\partial c}\mathbb{E}[\rho_{ \theta}(Y-c)] = 0 \Longleftrightarrow F(c) = \theta \Rightarrow Q_{Y}(\theta) = c \] Si ahora en lugar de usar la función de densidad de \(Y\) usamos la función de distribución empírica de \(Y\) que esta dada por:
\[ F_n(x) = \frac{1}{n} \sum_{k=1}^{n} 1_{\{Y_k \leq x\}} \] dónde \(\{Y_{i}\}_{i\geq 1}\) es una sucesión de v.a. tal que las v.a. \(Y_{i}\) se distribuye igual que \(Y\) para todo \(i\geq 1\) entonces el problema de minimización se convierte en
\[ \min_{c} \frac{1}{n} \sum_{k=1}^{n} \rho_{\theta}(y_k - c) \]
aquí, \(y_i\) son las observaciones de \(Y_i\). El valor de \(c\) que minimiza esta función es el cuantil \(\theta\) muestral (ya que las observaciones de las v.a. \(Y_{i}\) ya fueron observadas).
Supuesto de linealidad
Podemos extender la idea del cuantil como un problema de minimización a la regresión cuantílica. Dado que el cuantil \(\tau\)-ésimo (\(Q_{Y}(\tau)\)) minimiza la esperanza de la función de pérdida antes mencionada, el cuantil \(\tau\)-ésimo condicional a \(X\), que denotaremos por \(Q_{Y}(\tau|X)\), minimiza la esperanza condicional respecto a \(X\) de la función de pérdida antes mencionada.
Aún más, si asumimos \[Q_{Y}(\tau|x) = x^{T}\beta(\tau)\] es posible estimar a \(\beta(\tau )\) mediante el problema de minimización \[\min_{\beta}\sum_{k=1}^{n}\rho_{\theta}(y_{k}-x_{k}^{T}\beta(\tau))\]
El problema general de programación lineal
La programación lineal es una herramienta sumamente flexible, ampliamente aplicada en distintos campos y con diversos fines, tales como asignar recursos, programar trabajadores y planes, planificar inversiones en carteras, formular estrategias militares y definir estrategias de mercadeo. Un tratamiento exhaustivo del tema está fuera de este texto. Sin embargo, se presentan las ideas básicas útiles para comprender el procedimiento de resolución de la regresión cuantílica mediante programación lineal.
La programación lineal es una rama de la matemática que se ocupa de la asignación eficiente de recursos limitados a actividades conocidas con el objetivo de alcanzar una meta deseada, como minimizar costos o maximizar beneficios.
Las variables: \[ x_i \geq 0 \quad i=1, \ldots, n, \] cuyos valores deben decidirse de manera óptima, se denominan variables de decisión.
En un programa lineal general, el objetivo es encontrar un vector \(\mathbf{x}^* \in \mathbb{R}_{+}^n\) que minimice o maximice, el valor de una función lineal llamada función objetivo, dada entre todos los vectores \(\mathbf{x} \in \mathbb{R}_{+}^n\) que satisfacen un sistema dado de ecuaciones e inecuaciones lineales llamado sistema de restricciones. El papel de la linealidad es, por tanto, doble:
La linealidad de algunos modelos puede justificarse sobre la base de las propiedades típicas del problema. Sin embargo, algunos problemas no lineales pueden ser linealizados mediante un uso adecuado de transformaciones matemáticas. Este es el caso del problema de regresión cuantílica (QR, por sus siglas en inglés). La representación (o a veces la aproximación) de un problema mediante una formulación de programación lineal garantiza que existan procedimientos eficientes para calcular soluciones.
Para recapitular, un problema típico de programación lineal satisface las siguientes condiciones:
Geométricamente, esto restringe las soluciones a \(\mathbb{R}_{+}^n\). En caso de que el problema esté caracterizado por variables sin restricción de signo, es decir, variables que pueden ser positivas, negativas o cero, se puede utilizar un truco simple para restringirlas a variables no negativas, introduciendo dos variables no negativas \([x]^{+}\) y \([-x]^{+}\) y convirtiendo cualquier variable de decisión no restringida como la diferencia entre dos variables no negativas sin cambiar el problema de optimización: \[ \begin{aligned} & x=[x]^{+}-[-x]^{+} \\ & {[x]^{+} \geq 0} \\ & {[-x]^{+} \geq 0 .} \end{aligned} \]
Para \(x>0\), se fija \([x]^{+}=x\) y \([-x]^{+}=0\), mientras que para \(x<0\) se fija \([-x]^{+}=-x\) y \([x]^{+}=0\). Si \(x=0\), entonces \([x]^{+}=[-x]^{+}=0\).
La conversión de un problema de minimización a uno de maximización es trivial: maximizar \(z\) es equivalente a minimizar \(-z\).
Se pueden pasar las restricciones de una forma a otra. Por ejemplo, una restricción de desigualdad: \[ a_1 x_1+\cdots+a_i x_i+\cdots+a_n x_n \leq b \] puede convertirse en una restricción de mayor o igual simplemente multiplicándola por \(-1\). En cambio, puede convertirse en una restricción de igualdad añadiendo una variable no negativa (variable de holgura), \[ a_1 x_1+\cdots+a_i x_i+\cdots+a_n x_n+w=b, \quad w \geq 0 . \] Por otro lado, una restricción de igualdad, \[ a_1 x_1+\cdots+a_i x_i+\cdots+a_n x_n=b \] puede convertirse en una forma de desigualdad mediante la introducción de dos restricciones de desigualdad: \[ \begin{aligned} & a_1 x_1+\cdots+a_i x_i+\cdots+a_n x_n \leq b \\ & a_1 x_1+\cdots+a_i x_i+\cdots+a_n x_n \geq b . \end{aligned} \]
Esto significa que no hay diferencias en la forma en que se plantean las restricciones. Finalmente, vale la pena recordar que, desde un punto de vista geométrico, una ecuación lineal corresponde a un hiperplano, mientras que una desigualdad divide el espacio \(n\)-dimensional en dos, una parte en la que se satisface la desigualdad y otra en la que no.
Al combinar todas las desigualdades, el conjunto de vectores que satisfacen las restricciones, conjunto factible, es entonces la intersección de todos los semiespacios involucrados, los \(n\) semiespacios impuestos por la no negatividad de las variables de decisión y los \(m\) semiespacios impuestos por las \(m\) restricciones que regulan el proceso.
Consideremos un problema de programación lineal con \(n\) incógnitas (variables de decisión) y \(m\) restricciones. La formulación matricial del problema es la siguiente (forma estándar),
\[\begin{equation}\label{eq21} \begin{array}{ll} \operatorname{minimize} & \mathbf{c x} \\ \text { subject to } & \mathbf{A x} \leq \mathbf{b} \\ & \mathbf{x} \geq \mathbf{0} . \end{array} \end{equation}\]La forma estándar plantea las desigualdades como de tipo menor o igual y requiere la no negatividad de las variables de decisión.
Ya se ha mencionado cómo la condición \(\mathbf{x} \geq \mathbf{0}\) restringe a \(\mathbf{x}\) a \(\mathbb{R}_{+}^n\), es decir, al cuadrante positivo en el espacio \(n\)-dimensional. En \(\mathbb{R}^2\) es un cuarto del plano, en \(\mathbb{R}^3\) es un octavo del espacio, y así sucesivamente. Las restricciones \(\mathbf{A x} \leq \mathbf{b}\) generan \(m\) semiespacios adicionales. El conjunto factible consiste en la intersección de los \(m+n\) semiespacios mencionados. Dicho conjunto puede ser acotado, no acotado o vacío. La función de costo \(\mathbf{c x}\) genera una familia de hiperplanos paralelos, es decir, el hiperplano \(\mathbf{c x} = \text{constante}\) corresponde al conjunto de puntos cuyo costo es igual a dicha constante; cuando la constante varía, el hiperplano recorre todo el espacio \(n\)-dimensional. El óptimo \(\mathbf{x}^*\) es el vector \(n\)-dimensional, que garantiza el menor costo dentro del conjunto factible. Un problema que no tiene solución factible se llama no factible. Un problema con valores arbitrariamente grandes de la función objetivo se llama no acotado.
Asociado con cada programación lineal está su programa dual. Partiendo del problema lineal introducido anteriormente, también llamado el problema primal, su formulación dual comienza con la misma \(\mathbf{A}\) y \(\mathbf{b}\) pero invierte todo lo demás. El vector de costo primal \(\mathbf{c}\) y el vector de restricciones \(\mathbf{b}\) se intercambian con el vector de ganancia dual \(\mathbf{b}\) y \(\mathbf{c}\) se usa como el vector de restricciones. El desconocido dual (variable de decisión) \(\mathbf{y}\) es ahora un vector con \(m\) componentes, siendo las \(n\) restricciones representadas por \(\mathbf{A^{\top}y} \geq \mathbf{c}\). Dado el problema de programación lineal primal, el programa lineal dual asociado es, \[ \begin{array}{ll} \text{maximizar} & \mathbf{b^{\top} y} \\ \text{sujeto a} & \mathbf{A^{\top}y} \geq \mathbf{c} \\ & \mathbf{y} \geq \mathbf{0} . \end{array} \]
Los programas lineales vienen en pares primal/dual. Resulta que cada solución factible para uno de estos dos programas lineales da un límite sobre el valor objetivo óptimo para el otro.
La formulación de programación lineal para el problema QR
Para poder abordar el problema de la regresión cuantilica como un problema de programación lineal, se necesita llevar la función objetivo y el conjunto de restricciónes a una forma lineal. Para esto observe que el problema general se puede plantear como, \[ \begin{aligned} \min _{\boldsymbol{\beta}(\theta)}& \sum_{i=1}^n \rho_\theta\left(y_i-\mathbf{x}_i^{\top} \boldsymbol{\beta}(\theta)\right) \\ s.a&\quad \boldsymbol{y}=X\boldsymbol{\beta}(\theta)+\boldsymbol{\varepsilon}. \end{aligned} \] Para llevar este problema a uno que se pueda resolver en programación lineal, observe que la función objetivo está dada por, \[\rho_\theta\left(y_i-\boldsymbol{x}_i^{\top} \boldsymbol{\beta}(\theta)\right)= \begin{cases}\theta\left|y_i-\boldsymbol{x}_i^{\top} \boldsymbol{\beta}(\theta)\right| & \text { si } y_i-\boldsymbol{x}_i^{\top} \boldsymbol{\beta}(\theta)>0 \\ (1-\theta)\left|y_i-\boldsymbol{x}_i^{\top} \boldsymbol{\beta}(\theta)\right| & \text { si } y_i-\boldsymbol{x}_i^{\top} \boldsymbol{\beta}(\theta) \leq 0.\end{cases}\] Considere las siguientes funciones funciones, \[[y]_{+}=\left\{\begin{array}{ll}y & \text { si } y>0 \\ 0 & \text { o.c. }\end{array}, \quad[y]_{-}\left\{\begin{array}{lll}-y & \text { si } & y<0 \\ 0 & o . c\end{array}\right.\right.\] Note que si se utilizan estás funciones, la función objetivo se puede reescribir de la siguiente manera, \[\rho_\theta\left(y_i-\boldsymbol{x}_i^{\top} \boldsymbol{\beta}(\theta)\right)=\theta\left[y_i-\boldsymbol{x}_i^{\top} \boldsymbol{\beta}(\theta)\right]_{+}+(1-\theta)\left[\boldsymbol{x}_i^{\top} \boldsymbol{\beta}(\theta)-y_i\right]_{+}.\] Por otro lado, defina al vector de errores positivos \(\mathbf{s}_1\) y al vector de errores negativos \(\mathbf{s}_2,\) \[ \mathbf{s}_1=[\mathbf{y}-\mathbf{X} \boldsymbol{\beta}(\theta)]_{+}, \quad \mathbf{s}_2=[\mathbf{X} \boldsymbol{\beta}(\theta)-\mathbf{y}]_{+}. \] Note que \(\mathbf{s}_1\) y \(\mathbf{s}_2\) son dos vectores de dimensión n, donde \(\mathbf{s}_1\) en sus entradas solo tiene los errores positivos y donde los errores son negativos tiene cero, mientras que \(\mathbf{s}_2\) en sus entradas tiene el valor absoluto de los errores negativos y donde son negativos tienen cero. Observe que el problema original se puede reescribir de la siguiente manera, \[\begin{aligned} & \min _{\beta(\theta)} \sum_{i=1}^n \theta s_{1 i}+(1-\theta) s_{2 i} \\ \text { s.a } \boldsymbol{y} &=X \boldsymbol{\beta}(\theta)+\boldsymbol{s}_1-\boldsymbol{s}_2 \\ &=X [\boldsymbol{\beta}(\theta)]_{+}-X [\boldsymbol{\beta}(\theta)]_{-}+I \boldsymbol{s}_1-I \boldsymbol{s}_2.\end{aligned}\] Si se define a, \[ \boldsymbol{\psi} = \begin{pmatrix} [\boldsymbol{\beta}(\theta)]_{+} \\ [\boldsymbol{\beta}(\theta)]_{-} \\ \mathbf{s}_1 \\ \mathbf{s}_2 \end{pmatrix}, \quad \mathbf{d} = \begin{pmatrix} \mathbf{0}_p \\ \mathbf{0}_p \\ \theta \mathbf{1}_n \\ (1-\theta) \mathbf{1}_n \end{pmatrix}, \quad \mathbf{B} = \bigl[\,\mathbf{X}, -\mathbf{X}, \mathbf{I} , -\mathbf{I}\,\bigr], \] el problema primal se puede reescribir como, \[ \begin{aligned} \min & _\psi\quad \mathbf{d}^{\top} \psi \\ \text { sujeto a }& \quad\mathbf{B} \psi=\mathbf{y}. \end{aligned} \] De la expresión primal anterior se tiene que el problema dual está dado por, \[\begin{array}{ll}\max_{\boldsymbol{z}} & \boldsymbol{y}^{\top} \boldsymbol{z} \\ \text { s.a } & B^{\top} \boldsymbol{z} \leq 0 .\end{array}\] Observe que si se analiza el producto \(B^{\top} \boldsymbol{z},\) \[B^{\top} \boldsymbol{z}=\left[\begin{array}{c}X^{\top} \boldsymbol{z} \\ -X^{\top} \boldsymbol{z} \\ I \boldsymbol{z} \\ -I \boldsymbol{z}\end{array}\right] \leq\left[\begin{array}{c}\boldsymbol{0}_p \\ \boldsymbol{0}_p \\ \theta \boldsymbol{1}_n \\ (1-\theta) \boldsymbol{1}_n\end{array}\right] \Rightarrow \begin{aligned} & X^{\top} \boldsymbol{z}=0 \\ & \boldsymbol{z} \in[\theta-1, \theta]^n.\end{aligned}\] Esto debido a que en las dos primeras restricciones se tiene, \[ \begin{aligned} X^{\top} \boldsymbol{z} &\leq \boldsymbol{0}_p \\ X^{\top} \boldsymbol{z} &\geq \boldsymbol{0}_p \end{aligned} \] lo que equivale a la restricción de igualdad. Mientras que las dos ultimas restricciones,
\[ \begin{aligned} I \boldsymbol{z} &\leq \theta \boldsymbol{1}_n \\ I \boldsymbol{z} &\geq (\theta-1) \boldsymbol{1}_n, \end{aligned} \] lo que nos indica que cada entrada del vector \(\boldsymbol{z}\) debe ser menor o igual a \(\theta\) y que cada entrada del vector \(\boldsymbol{z}\) debe ser mayor o igual a \((\theta-1).\) Sea \(\boldsymbol{z}_0 \in [0,1]^n,\) observe que, \[\boldsymbol{z}_0 \in[0,1]^n \Rightarrow \boldsymbol{z}=\boldsymbol{z}_0+(\theta-1) \boldsymbol{1}_n.\] La función objetivo toma la forma, \[\max _{\boldsymbol{z}}\quad \boldsymbol{y}^{\top} \boldsymbol{z}=\max _{\boldsymbol{z}_0} y^{\top}\left(\boldsymbol{z}_0+(\theta-1) \boldsymbol{1}_{n}\right)=\max_{\boldsymbol{z}_0} \boldsymbol{y}^{\top} \boldsymbol{z}_0.\] Las restricciones se reescriben como, \[0=X^{\top} \boldsymbol{z}=X^{\top}\left(\boldsymbol{z}_0+(\theta-1) \boldsymbol{1}_n\right) \Leftrightarrow X^{\top} \boldsymbol{z}_0=(1-\theta) X^{\top} \boldsymbol{1}_n.\] Por lo tanto se tiene la siguiente formulación dual para el problema de regresión cuantílica, \[ \max _{\mathbf{z}}\left\{\mathbf{y}^{\top} \mathbf{z} \mid \mathbf{X}^{\top} \mathbf{z}=(1-\theta) \mathbf{X}^{\top} \mathbf{1}, \mathbf{z} \in[0,1]^n\right\}. \]
Principios de Equivarianza Básicos Existen algunas características que nos gustaría que preservara nuestro modelo, por ejemplo, si cambiáramos de escala el vector de respuesta, esperaríamos que las betas estimadas del la nueva regresión cuantílica se escalaran de la misma forma. De hecho esta propiedad se cumple y es una de las propiedades que enunciamos en el siguiente teorema. Antes que nada, denotamos a la estimación de \(\beta\) para el \(\tau\)-ésimo cuantil basado en las observaciones \((y,X)\) como \(\hat{\beta}(\tau; y,X)\).
Teorema: Sea \(A\) una matriz de \(p\times p\) no singular, \(\gamma\in \mathbb{R}^{n}\) y \(a>0\). Entonces, para cualquier \(\tau\in [0,1]\),
\(i)\quad \hat{\beta}(\tau;ay,X) = a\hat{\beta}(\tau;y, X)\)
\(ii) \quad \hat{\beta}(\tau, -ay,X) = -a\hat{\beta}(1-\tau; y,X)\)
\(iii) \quad\hat{\beta}(\tau;y+X\gamma,X) = \hat{\beta}(\tau;y, X)+\gamma\)
\(iv)\quad \hat{\beta}(\tau; y, XA) = A^{-1}\hat{\beta}(\tau:y, X)\)
Por poner un ejemplo, lo que estaría diciendo la primera probpiedad es: de hacer un escalamiento de la variable de respuesta, la estimación de la beta para el cuantíl \(\tau\) usando las covariables originales, sería lo mismo que multiplicar el vector estimado de las betas usando como variable de respuesta el vector original y las mismas covariables.
Equivarianza ante Transformaciones Monótonas
Sabemos que uno de los métodos más populares que se suelen ocupar para normalizar los datos de respuesta es aplicar la transformación de BoxCox, donde en muchos de los casos se toma la transformación \(Y_{i}^{\lambda} = log(Y_{i})\), uno de los casos más típicos es cuándo nuestros datos están sesgados hacia la izquierda y todos son positivos.Así, hacemos regresión sobre modelo
\[Y_{i}^{(\lambda)} = \beta_{0} + \beta_{1}x_{i,1}+\dots + \beta_{n}x_{i,n} + \varepsilon_{i} \]
así, asumimos que \(\mathbb{E}[Y^{(\lambda)}] = X^{T}\beta\) pero ocurre que no es posible recuperar el modelo original, en el sentido de que
\[ \exp\left[\mathbb{E}[Y_{i}^{(\lambda)}] \right] < \mathbb{E}[\exp\{\log(Y_{i})\}] = \mathbb{E}[Y_{i}]\] esto por la desigualdad de Jensen. Sin embargo, hay una propiedad de los cuantiles que sí nos permiten recuperar el modelo original para este caso y se enuncia a continuación.
Teorema: Sea \(h:\mathbb{R}\rightarrow \mathbb{R}\) una función monótona no decreciente en \(\mathbb{R}\), entonces, para cualquier variable aleatoria \(Y\) se tiene que \[Q_{h(Y)}(\tau) = h(Q_{Y}(\tau)).\] Dem. Este resultado se sigue facilmente de notar que
\[\tau = \mathbb{P}[Y\leq y_{\tau}] = \mathbb{P}[h(Y))\leq h(y_{\tau})]\] la última igualdad se da dado que \(h\) es una función monótona no decreciente. Así, por la últimaexpresión, se tiene que \[Q_{h(Y)}(\tau) = h(y_{\tau}) = h(Q_{Y}(\tau)).\]
Datos Censurados y Regresión Cuantílica
Una aplicación directa del teorema anterior y que además es una fortaleza de la regresión cuantílica es la aplicación a datos censurados. Considerese que se quiere hacer inferencia sobre el siguiente modelo \[y_{i} = x_{i}^{T}\beta + \varepsilon_{i}\] y supongamos que \(y_{i}\) no es directamente observable, que en su lugar observamos \(y^{*}_{i} = \max\{0,y_{i}\}\) (Censura por la izquierda). Entonces, del resultado anterior se sigue que \[Q_{y_{i}^{*}}(\tau|x) = \max\{0,x_{i}^{T}\beta(\tau)\}\] por lo que podemos encontrar a \(\beta_{0}\) y \(\beta_{1}\) mediante la solución de \[\min_{\beta}\sum_{k=1}^{n}\rho(y_{i}-\max \{0, x_{i}^{T}\beta\})\] También es posible hacer este análisis usando la función de máxima verosimilitud que nos llevaría a obtener los mismos estimadores que en la minimización del error cuadrático medio si la distirbución de los errores es Normal. En caso de que los erroes no sean normales, no se obtienen los mismos estimadores que minimizan el error cuadrático medio. Así, la virtud de la regresión cuantílica es que esta no asume una función de distirbución para los errores.
Ahora analizaremos el por qué bajo el criterio de la función de influencia de Hampel, la mediana es más robusta que la media. La función de influencia, introducida por Hampel, ofrece una descripción concisa de cómo un estimador \(\hat{\theta}\) evaluado en una distribución \(F\) se ve afectado al ``contaminar’’ \(F\). Formalmente, podemos ver a \(\hat{\theta}\) como una función de \(F\) y escribirlo como \(\hat{\theta}(F)\), y podemos considerar la contaminación de \(F\) reemplazando una pequeña cantidad de masa \(\varepsilon\) de \(F\) por una equivalente masa concentrada en \(y\), permitiéndonos escribir la función de distribución contaminada como
\[F_{\varepsilon}(x)=\varepsilon\delta_{y}(x)+(1-\varepsilon)F(x),\]
donde \(\delta_{y}\) denota la función de distribución que asigna masa 1 al punto \(y\). Expresamos la función de influencia de \(\hat{\theta}\) en \(F\) como
\[IF_{\hat{\theta}}(y,F)=\lim_{\varepsilon\to 0}\frac{\hat{\theta}(F_{ \varepsilon})-\hat{\theta}(F)}{\varepsilon}.\] Calculando la función de influencia de Hampel para la mediana se tiene que
\[\hat{\theta}(F_{\varepsilon})=\int x dF_{\varepsilon}(x)=\varepsilon y+(1- \varepsilon)\hat{\theta}(F)\]
y por lo tanto
\[IF_{\hat{\theta}}(y,F)=y-\hat{\theta}(F),\]
mientras que para la mediana,
\[\hat{\theta}(F_{\varepsilon})=F_{\varepsilon}^{-1}(1/2)\]
y es posible demostrar que
\[IF_{\hat{\theta}}(y,F)=\text{sgn}(y-\bar{\theta}(F))/f(F^{-1}(1/2)),\] De lo anterior sacamos las siguiente conclusiones
En el caso de la media, la influencia de contaminar \(F\) en \(y\) es proporcional a \(y\), lo que implica que una pequeña contaminación, , en un punto \(y\) suficientemente lejano a \(\theta(F)\) puede alejar arbitrariamente la media de su valor inicial en \(F\). Esto lo podemos ver dado que para una \(\varepsilon\) pequeña se tiene la aproximación \[\hat{\theta}(F_{\varepsilon}) - \hat{\theta}(F) \approx \epsilon (y-\hat{\theta}(F))\]
En contraste, la influencia de la contaminación en \(y\) sobre la mediana está por la constante \(1/f(F^{-1}(1/2))\), esto último dado que,
\[-\frac{\varepsilon}{f(F^{-1}(\frac{1}{2}))}\leq \hat{\theta}(F_{\varepsilon}) - \hat{\theta}(F) \leq \frac{\varepsilon}{f(F^{-1}(\frac{1}{2}))}\] aproximadamente.
Es posible extende esta misma idea a la regresión cuantílica, sin embargo, es posible demostrar que la función de información para la beta estimada en la regresión cuantílica no es sencible a lo extremos dela variable de respuesta, sin embargo, sí lo es para las covariables.
A continuación se introducirá la condición básica de optimalidad que caracteriza los cuantiles de regresión. Se ha visto que podemos restringir nuestra atención a soluciones candidatas de la forma básica, \[ b(h)=X(h)^{-1} y(h) \]
Para algunos \(h\), \(X(h)\) puede ser singular. Esto no debe preocuparnos, podemos restringirnos a \(b(h)\) con \(h \in \mathcal{H}^{*}=\{h \in \mathcal{H}:|X(h)| \neq 0\}\). También hemos visto que nuestra condición de optimalidad implica verificar que las derivadas direccionales sean no negativas en todas direcciones. Para verificar esto en \(b(h)\), debemos considerar, \[ \nabla R(b(h), w)=-\sum_{i=1}^{n} \psi_{\tau}^{*}\left(y_{i}-x_{i}^{\prime} b(h),-x_{i}^{\prime} w\right) x_{i}^{\prime} w, \] con \[ \psi^{*}(u, v)= \begin{cases}\tau-I(u<0) & \text{ si } u \neq 0 \\ \tau-I(v<0) & \text{ si } u=0\end{cases}. \] Reparametrizando las direcciones para que \(v=X(h) w\), tenemos optimalidad si y solo si,
\[ 0 \leq-\sum_{i=1}^{n} \psi_{\tau}^{*}\left(y_{i}-x_{i}^{\prime} b(h),-x_{i} X(h)^{-1} v\right) x_{i}^{\prime} X(h)^{-1} v \]
para todo \(v \in \mathbb{R}^{p}\). Notemos que para \(i \in h\), tenemos \(e_{i}^{\prime}=x_{i}^{\prime} X(h)^{-1}\), el \(i\)-ésimo vector unitario base de \(\mathbb{R}^{p}\), por lo que podemos reescribir esto como:
\[ 0 \leq-\sum_{i \in h} \psi_{\tau}^{*}\left(0, v_{i}\right) v_{i}-\xi^{\prime} v=-\sum_{i \in h}\left(\tau-I\left(v_{i}<0\right)\right) v_{i}-\xi^{\prime} v \] donde, \[ \xi\left(v_{i}\right)=\sum_{i \in \bar{h}} \psi_{\tau}^{*}\left(y_{i}-x_{i} b(h),-x_{i} X(h)^{-1} v_{i}\right) x_{i}^{\prime} X(h)^{-1} \]
Finalmente, notemos que el espacio de “direcciones” \(v \in \mathbb{R}^{p}\) está generado por \(v= \pm e_{k}, k=1, \ldots, p\). Es decir, la condición de derivada direccional se cumple para todo \(v \in \mathbb{R}^{p}\) si y solo si se cumple para las \(2p\) direcciones canónicas \(\{ \pm e_{i}: i=1, \ldots, p\}\). Así, para \(v=e_{i}\) tenemos las \(p\)-desigualdades:
\[ 0<-(\tau-1)+\xi_{i}\left(e_{i}\right) \quad i=1, \ldots, p \]
mientras que para \(v=-e_{i}\) tenemos:
\[ 0<\tau-\xi_{i}\left(-e_{i}\right) \quad i=1, \ldots, p \]
Combinando estas desigualdades, tenemos nuestra condición de optimalidad en toda su generalidad. Si ninguno de los residuos de las observaciones no básicas (\(i \in \bar{h}\)) es cero (como sería el caso con probabilidad uno si las \(y\) tienen una densidad respecto a la medida de Lebesgue), entonces la dependencia de \(\xi\) en \(v\) desaparece y podemos combinar los dos conjuntos de desigualdades para obtener:
\[ (\tau-1) 1_{p} \leq \xi_{h} \leq \tau 1_{p} \] Resumiendo la discusión anterior, podemos reformular el Teorema de Koenker y Bassett (1978) con ayuda de la siguiente definición introducida por Rousseeuw y Leroy (1987).
DEFINICIÓN: Diremos que las observaciones de regresión \((y, X)\) están en posición general si cualquier subconjunto de \(p\) de ellas produce un ajuste exacto único, es decir, para todo \(h \in \mathcal{H}^{*}\),
\[ y_{i}-x_{i} b(h) \neq 0 \quad \text{para todo } i \notin h. \]
Nótese que si las \(y_{i}\) tienen densidad respecto a la medida de Lebesgue, entonces las observaciones \((y, X)\) estarán en posición general con probabilidad uno.
TEOREMA. Si \((y, X)\) están en posición general, entonces existe una solución al problema de regresión cuantílica de la forma \(b(h)=X(h)^{-1} y(h)\) si y solo si para algún \(h \in \mathcal{H}^{*}\) se cumple
\[ (\tau-1) 1_{p} \leq \xi_{h} \leq \tau 1_{p} \tag{2.3.1} \]
donde \(\xi_{h}=\sum_{i \in \bar{h}} \psi_{\tau}\left(y_{i}-x_{i}^{\prime} b(h)\right) x_{i}^{\prime} X(h)^{-1}\) y \(\psi_{\tau}=\tau-I(u<0)\). Además, \(b(h)\) es la solución única si y solo si las desigualdades son estrictas; de lo contrario, el conjunto solución es la envolvente convexa de varias soluciones de la forma \(b(h)\).
Las estimaciones de QR pueden utilizarse para obtener intervalos de predicción de la distribución condicional de la variable respuesta. Este enfoque aprovecha la naturaleza no paramétrica de QR y ofrece una herramienta flexible para estimar intervalos de predicción en cualquier valor específico de las variables explicativas. A diferencia de los intervalos de predicción ordinarios utilizados en regresión por mínimos cuadrados, estos intervalos no son sensibles a desviaciones de los supuestos distribucionales del modelo lineal normal clásico.
Para un modelo dado, el intervalo proporcionado por dos estimaciones cuantílicas distintas, \(\hat{q}_Y(\theta_1, X=x)\) y \(\hat{q}_Y(\theta_2, X=x)\), en cualquier valor específico del regresor \(X\), constituye un intervalo de predicción del \((\theta_2-\theta_1)\%\) para una única observación futura. Si consideramos, por ejemplo, \(\theta_1=0.1\) y \(\theta_2=0.9\), entonces el intervalo \([\hat{q}_Y(\theta_1=0.1, X=x); \hat{q}_Y(\theta_2=0.9, X=x)]\) puede interpretarse como un intervalo de predicción del 80% correspondiente al valor dado \(X=x\).
El uso de este tipo de intervalo de predicción garantiza un nivel de cobertura consistente con el nivel nominal independientemente de la naturaleza del término de error, ya que QR no se ve afectado por desviaciones de los supuestos clásicos del modelo lineal normal. Además, este nivel de cobertura se mantiene estable para cualquier valor específico de \(x\). Esta es una gran ventaja respecto a la regresión clásica, donde la correspondencia entre los niveles de cobertura empíricos y nominales solo se garantiza en el caso de los supuestos distribucionales usuales (modelo \(_1\)), mientras que se deteriora cuando aumenta el grado de violación de estos supuestos.
Supongamos que \(Y_{1}, \ldots, Y_{n}\) son variables aleatorias independientes e idénticamente distribuidas (iid) con función de distribución común \(F\) y asumamos que \(F\) tiene una densidad continua \(f\) en una vecindad de \(\xi_{\tau}=F^{-1}(\tau)\) con \(f\left(\xi_{\tau}\right)>0\). La función objetivo del \(\tau\)-ésimo cuantil muestral
\[ \hat{\xi}_{\tau} \equiv \inf _{\xi}\left\{\xi \in \mathbb{R} \mid \sum \rho_{\tau}\left(Y_{i}-\xi\right)=\min \right\} \]
Observe que, \[\rho_{\tau}(y) = [\tau-I(y<0)]|y| = |y|(1-\theta)I(y\leq 0) + |y|\theta I(y>0),\]
luego la función objetivo se puede ver como la suma de funciones convexas y por lo tanto es convexa. Por otro lado, recuerde que, \[\begin{aligned} & \rho_\tau\left(Y_i-\xi\right)= \begin{cases}\tau\left(Y_i-\xi\right), & \text { si }\left(Y_i-\xi\right)>0 \\ (1-\tau)(-1)\left(Y_i-\xi\right) & \text { si }\left(Y_{i}- \xi\right)<0\end{cases} \\ & \rho_\tau^{\prime}\left(Y_i-\xi\right)= \begin{cases}-\tau & \text { si }\left(Y_i-\xi\right)>0 \\ (1-\tau) & \text { si }\left(Y_i-\xi\right)<0.\end{cases}, \end{aligned}\] de este modo, la derivada respecto a \(\xi\) se puede escribir como, \[\frac{\partial\rho_\tau^\prime(Y_i-\xi)}{\partial \xi}=-\tau I_{Y_i-\xi>0}+(1-\tau)I_{Y_i-\xi<0}=I_{Y_i<\xi}-\tau \] Se sigue, que el gradiente tiene la forma,
\[ g_{n}(\xi)=\sum_{i=1}^{n}\left(I\left(Y_{i}<\xi\right)-\tau\right) \] Luego, por la monotonía del gradiente en funciones convexas se tiene,
\[ \begin{aligned} P\left\{\hat{\xi}_{\tau}>\xi\right\} & =P\left\{g_{n}(\xi)<0\right\} \\ & =P\left\{\sum I\left(Y_{i}<\xi\right)<n \tau\right\} \\ & =P\{B(n, F(\xi))<n \tau\}, \end{aligned} \]
donde \(B(n, p)\) denota una variable aleatoria binomial con parámetros \((n, p)\). Así, dejando que \(m=\lceil n \tau\rceil\) denote el entero más pequeño \(\geq n \tau\), podemos expresar la función de distribución \(G(\xi) \equiv P\left\{\hat{\xi}_{\tau} \leq \xi\right\}\) de \(\hat{\xi}(\tau)\) usando la función beta incompleta como
\[ \begin{aligned} G(\xi) & =\sum_{k=m}^{n}\binom{n}{k} F(\xi)^{k}(1-F(\xi))^{n-k} \\ & =n\binom{n-1}{m-1} \int_{0}^{F(\xi)} t^{m-1}(1-t)^{n-m} d t \end{aligned} \]
Diferenciando, obtenemos la función de densidad para \(\hat{\xi}(\tau)\)
\[ g(\xi)=n\binom{n-1}{m-1} F(\xi)^{m-1}(1-F(\xi))^{n-m} f(\xi) \]
Esta forma de la densidad puede deducirse directamente notando que el evento \(\{x<\) \(\left.Y_{(m)}<x+\delta\right\}\) requiere que \(m-1\) observaciones estén por debajo de \(x\), \(n-m\) estén por encima de \(x+\delta\) y una esté en el intervalo \((x, x+\delta)\). El número de formas en que puede ocurrir esta disposición es
\[ \frac{n!}{(m-1)!1!(n-m)!}=n\binom{n-1}{m-1} \]
y cada disposición tiene la probabilidad \(\left.F(\xi)^{m-1}(1-F(\xi))^{n-m}[F(\xi+\delta)-F(\xi))\right]\). Así,
\[ P\left\{x<Y_{(m)}<x+\delta\right\}=n\binom{n-1}{m-1} F(\xi)^{m-1}(1-F(\xi))^{n-m} f(\xi) \delta+o\left(\delta^{2}\right), \]
y obtenemos la función de densidad dividiendo ambos lados por \(\delta\) y dejando que tienda a cero.\ Este enfoque también puede usarse para construir intervalos de confianza para \(\xi_{\tau}\) de la forma
\[ P\left\{\hat{\xi}_{\tau_{1}}<\xi_{\tau}<\hat{\xi}_{\tau_{2}}\right\}=1-\alpha \]
donde \(\tau_{1}\) y \(\tau_{2}\) se eligen para satisfacer
\[ P\left\{n \tau_{1}<B(n, \tau)<n \tau_{2}\right\}=1-\alpha \]
En el caso de \(F\) continua, estos intervalos tienen la notable propiedad de ser libres de distribución, es decir, se mantienen independientemente de \(F\). En el caso de \(F\) discreta, se pueden construir límites libres de distribución estrechamente relacionados para \(P\left\{\hat{\xi}_{\tau_{1}}<\xi_{\tau}<\hat{\xi}_{\tau_{2}}\right\}\) y \(P\left\{\hat{\xi}_{\tau_{1}} \leq \xi_{\tau} \leq \hat{\xi}_{\tau_{2}}\right\}\).
La densidad de muestra finita del estimador de regresión cuantílica \(\hat{\beta}(\tau)\) bajo errores iid puede derivarse de la condición de subgradiente introducida anteriormente de una manera muy similar a la derivación de la densidad en el caso de una muestra dado anteriormente.
Teorema de (Bassett y Koenker (1978)). Considere el modelo lineal
\[ Y_{i}=x_{i}^{\prime} \beta+u_{i} \quad i=1, \ldots, n, \]
con errores iid \(\left\{u_{i}\right\}\) que tienen función de distribución común \(F\) y densidad estrictamente positiva \(f\) en \(F^{-1}(\tau)\). Entonces la densidad de \(\hat{\beta}(\tau)\) toma la forma,
\[ g(b)=\sum_{h \in \mathcal{H}} P\left\{\xi_{h}(b) \in \mathcal{C}\right\}|X(h)| \prod_{i \in h} f\left(x_{i}^{\prime}(b-\beta(\tau))+F^{-1}(\tau)\right) \tag{3.1.2} \]
donde \(\xi_{h}(b)=\sum_{i \in \bar{h}} \psi_{\tau}\left(y_{i}-x_{i} b\right) x_{i}^{\prime} X(h)^{-1}\) y \(\mathcal{C}\) denota el cubo \([\tau-1, \tau]^{p}\).\ A continuación se da una demostración de este resultado. Bajo las condiciones de soluciones básicas se sabe que, \(\hat{\beta}(\tau)=b(h) \equiv(X(h))^{-1} Y(h)\) si y solo si \(\xi_{h}(b(h)) \in\) \(\mathcal{C}\). Para cualquier \(b \in \mathbb{R}^{p}\), sea \(B(b, \delta)=b+[-\delta / 2, \delta / 2]^{p}\) el cubo centrado en \(b\) con aristas de longitud \(\delta\) y escribimos
\[ \begin{align*} P\{\beta(\tau) \in B(b, \delta) & =\sum_{h \in \mathcal{H}} P\left\{b(h) \in B(b, \delta), \xi_{h}(b(h)) \in \mathcal{C}\right\} \tag{3.1.3}\\ & =\sum_{h \in \mathcal{H}} E I(b(h) \in B(b, \delta)) P\left\{\xi_{h}(b(h)) \in \mathcal{C} \mid Y(h)\right\}. \end{align*} \]
La probabilidad condicional anterior se define a partir de la distribución de \(\xi_{h}(b(h))\) que es una variable aleatoria discreta (tomando \(2^{n-p}\) valores para cada \(h \in \mathcal{H}\) ). Cuando \(\delta \rightarrow 0\), esta probabilidad condicional tiende a \(P\left\{\xi_{h}(b) \in \mathcal{C}\right\}\); esta probabilidad es independiente de \(Y(h)\).
Ahora, dividimos ambos lados por Volumen \((B(b, \delta))=\delta^{p}\) y dejamos \(\delta \rightarrow 0\) para expresarlo en forma de densidad. La probabilidad condicional tiende a \(P\left\{\xi_{h}(b) \in \mathcal{C}\right\}\) que ya no depende de \(Y(h)\). El otro factor tiende a la densidad conjunta del vector \(X(h)^{-1} Y(h)\) que como
\[ f_{Y(h)}(y)=\prod_{i \in h} f\left(y_{i}-x_{i}^{\prime} \beta\right) \]
puede escribirse como
\[ \begin{align*} f_{(X(h))^{-1} Y(h)}(b) & =|X(h)| \prod_{i \in h} f\left((X(h) b)_{i}-x_{i}^{\prime} \beta(\tau)+F^{-1}(\tau)\right) \tag{3.1.4}\\ & =|X(h)| \prod_{i \in h} f\left(x_{i}^{\prime}(b-\beta(\tau))+F^{-1}(\tau)\right) \end{align*} \]
Nótese que los \(h\) para los cuales \(X(h)\) es singular no contribuyen a la densidad. El resultado ahora sigue al reensamblar las piezas.
Desafortunadamente, desde un punto de vista práctico los \(\binom{n}{p}\) sumandos no son muy manejables en la mayoría de las aplicaciones y, como en la teoría de mínimos cuadrados, debemos recurrir a aproximaciones asintóticas para una teoría de distribución adaptable a la inferencia estadística práctica. En la siguiente sección tratamos de proporcionar una introducción heurística a la teoría asintótica de la regresión cuantílica.
La visión de optimización de los cuantiles muestrales también permite un enfoque elemental de su teoría asintótica. Supongamos que \(Y_{1}, \ldots, Y_{n}\) son independientes e idénticamente distribuidos (iid) de la distribución \(F\) y para algún cuantil \(\xi_{\tau}=F^{-1}(\tau)\) asumamos que \(F\) tiene una densidad continua \(f\) en \(\xi_{\tau}\) con \(f\left(\xi_{\tau}\right)>0\). La función objetivo del \(\tau\)-ésimo cuantil muestral
\[ \hat{\xi}_{\tau} \equiv \inf _{\xi}\left\{\xi \in \mathbb{R} \mid \sum \rho_{\tau}\left(Y_{i}-\xi\right)=\min !\right\} \]
como ya hemos notado, esta función objetivo es la suma de funciones convexas, por lo tanto es ella misma convexa. Consecuentemente su gradiente
\[ g_{n}(\xi)=n^{-1} \sum_{i=1}^{n}\left(I\left(Y_{i}<\xi\right)-\tau\right) \]
es monótono en \(\xi\). Por supuesto, cuando \(\xi\) iguala a uno de los \(Y_{i}\) entonces este “gradiente” necesita la interpretación de subgradiente discutida anteriormente, pero esto no es crucial para el argumento que sigue. Por monotonía, \(\hat{\xi}_{\tau}\) es mayor que \(\xi\) si y solo si \(g_{n}(\xi)<0\). Así que
\[ \begin{aligned} P\left\{\sqrt{n}\left(\hat{\xi}_{\tau}-\xi_{\tau}\right)>\delta\right\} & =P\left\{\left(\hat{\xi}_{\tau}\right)>\frac{\delta}{\sqrt{n}}+\xi_\tau\right\} \\&=P\left\{g_{n}\left(\xi_{\tau}+\delta / \sqrt{n}\right)<0\right\} \\ & =P\left\{n^{-1} \sum\left(I\left(Y_{i}<\xi_{\tau}+\delta / \sqrt{n}\right)-\tau\right)<0\right\} \end{aligned} \]
Así, hemos reducido el comportamiento de \(\hat{\xi}_{\tau}\) a un problema de teorema del límite central DeMoivre-Laplace en el que tenemos un arreglo de variables aleatorias Bernoulli. Los sumandos toman los valores \((1-\tau)\) y \(-\tau\) con probabilidades \(F\left(\xi_{\tau}+\delta / \sqrt{n}\right)\) y \(1-F\left(\xi_{\tau}+\delta / \sqrt{n}\right)\). Dado que
\[ E [g_{n}\left(\xi_{\tau}+\delta / \sqrt{n}\right)]=\left(F\left(\xi_{\tau}+\delta / \sqrt{n}\right)-\tau\right) \rightarrow f\left(\xi_{\tau}\right) \delta / \sqrt{n} \]
y
\[ \begin{aligned} &V\left(g_{n}\left(\xi_{\tau}+\delta / \sqrt{n}\right)\right)=\\ &F\left(\xi_{\tau}+\delta / \sqrt{n}\right)\left(1-F\left(\xi_{\tau}+\delta / \sqrt{n}\right)\right) / n \rightarrow \tau(1-\tau) / n . \end{aligned}\]
podemos establecer \(\omega^{2}=\tau(1-\tau) / f^{2}\left(\xi_{\tau}\right)\) y escribir
\[ \begin{aligned} P\left\{\sqrt{n}\left(\hat{\xi}_{\tau}-\xi_{\tau}\right)>\delta\right\} & =P\left\{\frac{g_{n}\left(\xi_{\tau}+\delta / \sqrt{n}\right)-f\left(\xi_{\tau}\right) \delta / \sqrt{n}}{\sqrt{\tau(1-\tau) / n}}<-\omega^{-1} \delta\right\} \\ & \rightarrow 1-\Phi\left(\omega^{-1} \delta\right) \end{aligned} \]
y por lo tanto
\[ \sqrt{n}\left(\hat{\xi}_{\tau}-\xi_{\tau}\right) \leadsto \mathcal{N}\left(0, \omega^{2}\right) \tag{3.2.1} \]
El efecto \(\tau(1-\tau)\) tiende a hacer \(\hat{\xi}_{\tau}\) más preciso en las colas, pero esto típicamente sería dominado por el efecto del término de densidad que tiende a hacer \(\hat{\xi}_{\tau}\) menos preciso en regiones de baja densidad.
Extendiendo el argumento anterior para considerar la forma límite de la distribución conjunta de varios cuantiles, sea \(\hat{\zeta}_{n}=\left(\hat{\xi}_{\tau_{1}}, \ldots, \hat{\xi}_{\tau_{m}}\right)\) con \(\zeta_{n}=\left(\xi_{\tau_{1}}, \ldots, \xi_{\tau_{m}}\right)\) y obtenemos
\[ \sqrt{n}\left(\hat{\zeta}_{n}-\zeta\right) \leadsto \mathcal{N}(0, \Omega) \tag{3.2.2} \]
donde \(\Omega=\left(\omega_{i j}\right)=\left(\tau_{i} \wedge \tau_{j}-\tau_{i} \tau_{j}\right) /\left(f\left(F^{-1}\left(\tau_{i}\right) f\left(F^{-1}\left(\tau_{j}\right)\right)\right.\right.\). Estos resultados para los cuantiles muestrales ordinarios en el modelo de una muestra se generalizan elegantemente al modelo de regresión lineal clásico
\[ y_{i}=x_{i}^{\prime} \beta+u_{i} \]
con errores iid \(\left\{u_{i}\right\}\). Supongamos que los \(\left\{u_{i}\right\}\) tienen función de distribución común \(F\) con densidad asociada \(f\), con \(f\left(F^{-1}\left(\tau_{i}\right)\right)>0\) para \(i=1, \ldots, m\), y \(n^{-1} \sum x_{i} x_{i}^{\prime} \equiv Q_{n}\) converge a una matriz definida positiva \(Q_{0}\). Entonces la distribución asintótica conjunta de los estimadores de regresión cuantílica \(m p\)-variantes \(\hat{\zeta}_{n}=\left(\hat{\beta}_{n}\left(\tau_{1}\right)^{\prime}, \ldots, \hat{\beta}_{n}\left(\tau_{m}\right)^{\prime}\right)^{\prime}\) toma la forma
\[ \begin{equation*} \sqrt{n}\left(\hat{\zeta}_{n}-\zeta\right)=\left(\sqrt{n}\left(\hat{\beta}_{n}\left(\tau_{j}\right)-\beta\left(\tau_{j}\right)\right)\right)_{j=1}^{m}=\mathcal{N}\left(0, \Omega \otimes Q_{0}^{-1}\right) . \tag{3.2.3} \end{equation*} \]
En el escenario de regresión con errores iid, la forma de los \(\beta\left(\tau_{j}\right)\) es particularmente simple como hemos visto. Los planos cuantiles condicionales de \(y \mid x\) son paralelos, así que suponiendo que la primera coordenada de \(\beta\) corresponde al parámetro “intercepto”, tenemos \(\beta(\tau)=\beta+\xi_{\tau} e_{1}\) donde \(\xi_{\tau}=F^{-1}(\tau)\) y \(e_{1}\) es el primer vector unitario base de \(\mathbb{R}^{p}\). Dado que \(\Omega\) toma la misma forma que en el escenario de una muestra, muchos de los resultados clásicos sobre L-estadísticas pueden llevarse directamente a la regresión con errores iid usando este resultado.
Interpretación de pendientes
Visualización progresiva
Curvas cuantílicas para \(\tau = 0.1, 0.25, 0.5, 0.75, 0.9\).
Simular datos
set.seed(42)
n <- 500
education <- sample(5:20, n, replace = TRUE)
income <- 5 + 2.5 * education + (education - 12) * rnorm(n, sd = 3)
df <- data.frame(education = education, income = income)
# Agregar constante (intercepto) para la regresión
X <- model.matrix(~ education, data = df)Función de visualización
library(quantreg)
library(ggplot2)
set.seed(42)
n <- 500
# Función para graficar regresión cuantílica para un valor dado de tau
plot_quantile_regression <- function(df, tau = 0.5) {
model <- rq(income ~ education, data = df, tau = tau)
df$pred <- predict(model)
ggplot(df, aes(x = education, y = income)) +
geom_point(alpha = 0.3, color = "gray40") +
geom_line(aes(y = pred), color = "red", size = 1) +
labs(
title = paste("Regresión Cuantílica para τ =", tau),
x = "Años de Educación",
y = "Ingreso estimado"
) +
theme_minimal()
}Interactivo
```{python}
from ipywidgets import interact, FloatSlider
def f(x):
return x
interact(f, x=10);
#(plot_quantile_regression, tau=FloatSlider(value=0.5, min=0.05, max=0.95, step=0.05))
```Distribución de ingresos: Estudio ENIGH 2022 (CDMX): ingreso vs. educación.
En el análisis realizado con los datos de la ENIGH 2022 para la Ciudad de México, se aplicó un modelo de regresión cuantíl para estimar el efecto de la educación del jefe del hogar sobre el ingreso trimestral. Los resultados muestran una fuerte heterogeneidad en los retornos educativos a lo largo de la distribución del ingreso. Específicamente, mientras que en el cuantil 5% el coeficiente asociado a un año adicional de educación es cercano a $2,000 MXN, en el cuantil 95% asciende a más de $19,000 MXN. Este patrón evidencia que el efecto de la educación no es constante, sino que se amplifica significativamente en los hogares de mayores ingresos.
Estos hallazgos coinciden con lo reportado en la literatura económica. Por ejemplo, Koenker y Hallock (2001) documentan que los retornos a la educación tienden a ser mayores en los cuantiles superiores del ingreso, lo que refleja la interacción entre capital humano y segmentación del mercado laboral. Asimismo, estudios aplicados al caso mexicano, como los de López-Acevedo (2001) y Campos-Vázquez (2013), confirman que la educación produce mayores rendimientos en ocupaciones mejor remuneradas, mientras que en los estratos bajos su efecto es limitado por restricciones estructurales como la informalidad, la baja productividad o el acceso desigual a oportunidades laborales.
En conjunto, estos resultados refuerzan la necesidad de adoptar políticas complementarias a la expansión educativa, tales como mejoras en la calidad del empleo, programas de capacitación focalizados y estrategias para reducir la segmentación del mercado laboral. La regresión cuantíl permite visualizar estas desigualdades de forma más precisa que los métodos tradicionales, al capturar no solo el efecto promedio, sino la forma en que dicho efecto varía en distintos contextos sociales y económicos.
Resultados Empíricos
Visualización interactiva con Plotly
```{python}
def plot_quantile_plotly(tau=0.5):
X = sm.add_constant(df['educacion_jefe'])
y = df['ingreso_trimestral']
model = sm.QuantReg(y, X).fit(q=tau)
beta = model.params['educacion_jefe']
pred = model.predict(X)
fig = go.Figure()
# Puntos
fig.add_trace(go.Scatter(
x=df['educacion_jefe'], y=y,
mode='markers',
name='Datos',
marker=dict(color='rgba(0,0,0,0.3)', size=5)
))
# Línea de regresión cuantíl
fig.add_trace(go.Scatter(
x=df['educacion_jefe'], y=pred,
mode='lines',
name=f'Regresión Cuantíl (τ = {tau:.2f})',
line=dict(color='red', width=3)
))
fig.update_layout(
title=f"CDMX - ENIGH 2022 - Regresión Cuantíl (τ = {tau:.2f})",
xaxis_title="Educación del jefe del hogar",
yaxis_title="Ingreso trimestral (MXN)",
template='plotly_white'
)
fig.add_annotation(
text=f"β = {beta:.2f}",
xref="paper", yref="paper",
x=0.05, y=0.95,
showarrow=False,
font=dict(size=14, color="red"),
bgcolor="rgba(255,255,255,0.7)",
bordercolor="red",
borderwidth=1
)
fig.show()
``````{python}
# --- Control deslizante interactivo ---
interact(plot_quantile_plotly, tau=FloatSlider(value=0.5, min=0.05, max=0.95, step=0.05))
```Coeficientes por cuantíl
| Cuantil \(\tau\) | \(\hat\beta_\tau\) (MXN/año) |
|---|---|
| 0.05 | 2 000 |
| 0.50 | 5 300 |
| 0.95 | 19 000 |
Interpretación
Política pública
Métodos de validación en regresión cuantíl
A diferencia de la regresión lineal ordinaria, donde se evalúa el ajuste del modelo mediante la suma de cuadrados de los residuos y el coeficiente de determinación \(R^2\), la regresión cuantíl utiliza criterios diferentes, ya que no modela la media condicional sino los cuantiles condicionales de la variable respuesta. Por ello, los métodos de validación deben adaptarse a la función de pérdida que optimiza este tipo de modelos: la pérdida de cuantíl o función de “check”.
Una primera forma de validación consiste en verificar que los residuos estén distribuidos conforme al cuantil modelado. Por ejemplo, si se estima el cuantil \(\tau = 0.25\), se espera que aproximadamente el 25% de los residuos \(e_i = y_i - \hat{Q}_\tau(x_i)\) sean negativos. Esta proporción sirve como un chequeo rápido de la coherencia del modelo, y puede interpretarse como una forma empírica de validar si la regresión cuantíl está ubicando correctamente la frontera de probabilidad deseada.
Un segundo criterio más formal es el pseudo \(R^1(\tau)\), propuesto por Koenker y Machado (1999). Este coeficiente se define como una analogía al \(R^2\) clásico, pero adaptado a la pérdida asimétrica de la regresión cuantíl:
\[ R^1(\tau) = 1 - \frac{\sum_{i=1}^n \rho_\tau(y_i - x_i^\top \hat{\beta}_\tau)}{\sum_{i=1}^n \rho_\tau(y_i - \tilde{y})} \]
donde \(\rho_\tau(u) = u(\tau - \mathbf{1}_{\{u < 0\}})\) es la función de pérdida de cuantíl, y \(\tilde{y}\) es el cuantil empírico \(\tau\) de la muestra, es decir, el valor que se obtendría al predecir siempre el mismo cuantíl sin usar ninguna covariable.
Este índice compara el desempeño del modelo completo con covariables contra un modelo base sin explicativas. Si el valor es cercano a 1, el modelo está explicando bien el cuantil condicional; si es cercano a 0, las covariables apenas aportan; y si es negativo, el modelo con variables explicativas resulta peor que simplemente usar el cuantil empírico de la muestra.
En conjunto, estas herramientas permiten validar modelos de regresión cuantíl respetando su estructura no paramétrica y robusta, y ofrecen criterios objetivos para comparar diferentes especificaciones de modelo o para identificar efectos heterogéneos mal capturados.
Residuos
\[ e_i^{(\tau)} = y_i - x_i^\top \hat\beta_\tau. \] Aproximadamente \(\tau\cdot100\%\) residuos negativos.
Pseudo \(R^1(\tau)\)
\[ R^1(\tau) = 1 - \frac{\sum_i \rho_\tau(y_i - x_i^\top \hat\beta_\tau)} {\sum_i \rho_\tau(y_i - \tilde y)}, \] donde \(\tilde y\) es el cuantil empírico de \(y\) sin covariables.
Diagnóstico gráfico
Gráfico Q-Q empírico para validación de residuos
Otra herramienta útil para diagnosticar la calidad del ajuste en regresión cuantíl es el gráfico Q-Q (cuantil-cuantil) empírico de los residuos. A diferencia del gráfico Q-Q tradicional en regresión lineal, que compara los residuos contra una distribución normal teórica, en regresión cuantíl se utiliza un enfoque no paramétrico y empírico: se comparan los residuos ordenados contra los cuantiles teóricos esperados de una distribución simétrica para el nivel de \(\tau\) considerado.
Recordemos que si el modelo está bien especificado, los residuos \(e_i = y_i - x_i^\top \hat{\beta}_\tau\) deberían estar distribuidos de forma tal que aproximadamente un \(\tau \cdot 100\%\) de ellos sean negativos, y el resto positivos. Esta distribución esperada no es necesariamente normal, pero sí debe mantener una simetría alrededor del cero para ciertos cuantiles centrales (por ejemplo, la mediana). Por eso, en lugar de usar una referencia gaussiana, se puede comparar la distribución de los residuos observados con la de una distribución simétrica simulada o simplemente contrastar su comportamiento con lo esperado bajo el supuesto de independencia e identidad del modelo.
En la práctica, el gráfico Q-Q empírico se construye de la siguiente manera:
Un buen ajuste se refleja en una nube de puntos alineada con la diagonal \(y = x\). Desviaciones sistemáticas indican asimetrías, colas más pesadas o estructuras no capturadas por el modelo. Este tipo de diagnóstico visual es especialmente útil cuando se combinan regresiones cuantílicas para distintos valores de \(\tau\), ya que permite comparar cómo cambia la forma de los residuos a lo largo de la distribución condicional.
¿Cómo se construye un Q-Q plot en regresión cuantíl?
En un Q-Q plot tradicional (como en regresión lineal), transformamos los datos usando la inversa de una distribución teórica —por ejemplo, la normal— y los comparamos con los cuantiles empíricos ordenados. Sin embargo, en regresión cuantíl este enfoque no aplica directamente porque:
Enfoque 1: Q-Q contra cuantiles simulados bajo el modelo
Se calculan los residuos:
\[
e_i = y_i - x_i^\top \hat{\beta}_\tau
\]
Se ordenan de menor a mayor:
\[
e_{(1)} < e_{(2)} < \cdots < e_{(n)}
\]
Se calculan cuantiles teóricos esperados bajo la hipótesis de un modelo bien ajustado: No hay una distribución explícita, pero se puede construir una referencia simétrica por remuestreo (bootstrap) o suponer una estructura empírica.
Se grafican los cuantiles teóricos vs. los residuos ordenados, como en un Q-Q clásico.
Este enfoque busca detectar si los residuos muestran asimetrías o estructuras no capturadas por el modelo, sin requerir una distribución paramétrica.
Enfoque 2: Q-Q de residuos contra distribución uniforme centrada
Otra estrategia consiste en analizar los residuos transformados como si siguieran una distribución uniforme si el modelo está bien especificado:
Calcular los residuos \(e_i\).
Asignarles sus posiciones empíricas (ranks normalizados):
\[
u_i = \frac{\text{rango}(e_i)}{n+1}
\]
Compararlos con los cuantiles de una distribución uniforme \(U(0,1)\).
Este enfoque es más informal pero útil para validar si la dispersión y el ordenamiento de los residuos parecen razonables.
Idea central
En regresión cuantíl, no comparamos contra una distribución normal teórica, sino que evaluamos si los residuos tienen la forma esperada dados los cuantiles estimados. Esto permite detectar asimetrías, colas pesadas o falta de ajuste sin imponer una estructura paramétrica estricta.
Referencias:
statsmodels en Python.Equipo 3:
Omar García
Rodrigo Gozaga
Joksan Avendaño
Modelo lineal simple
En el caso usual, el paramétrico, la regresión lineal simple según @olaya2020 es uno de los modelos más comunes. Un modelo estadístico lineal general se escribe como:
\[ y = f(x) + \varepsilon \tag{1}\label{eq:modelo} \]
donde \(y\) y \(\varepsilon\) son variables aleatorias, y \(f\) es una función determinística de una variable no aleatoria.
En un modelo lineal, \(f\) es lineal en los parámetros desconocidos:
\[ f(x) = \beta_0 + \beta_1 x \]
Así, el modelo se expresa como: \[ y = \beta_0 + \beta_1 x + \varepsilon, \quad \mathbb{E}[\varepsilon] = 0, \quad \text{Var}(\varepsilon) = \sigma^2 \]
Los estimadores de mínimos cuadrados para \(\beta_0\) y \(\beta_1\) son: \[ \hat{\beta}_1 = \frac{\sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y})}{\sum_{i=1}^n (x_i - \bar{x})^2}, \quad \hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \bar{x} \]
El modelo estimado es entonces:
\[ \hat{f}(x) = \hat{\beta}_0 + \hat{\beta}_1 x \]
Para los siguientes datos simulados, se calcula y grafica la regresión de los datos:
#>
#> Call:
#> lm(formula = y ~ x)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -4.6518 -2.0117 -0.7275 1.8043 9.9492
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -0.6421 0.6961 -0.922 0.359
#> x 2.0113 0.1142 17.609 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 2.997 on 98 degrees of freedom
#> Multiple R-squared: 0.7598, Adjusted R-squared: 0.7574
#> F-statistic: 310.1 on 1 and 98 DF, p-value: < 2.2e-16| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | -0.642 | 0.696 | -0.922 | 0.359 |
| x | 2.011 | 0.114 | 17.609 | 0.000 |
library(ggplot2)
ggplot(data.frame(x, y), aes(x, y)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE, color = "red") +
ggtitle("Ajuste lineal")
Como se puede observar en Figura 7.1 y en tab:summary-modelo-l el modelo para ser no muy adecuado, esto se discutira mas adelante.
Modelos lineal múltiple
Supongamos que tenemos datos \((x_1, Y_1), \ldots, (x_n, Y_n)\) donde \(Y_i \in \mathbb{R}\) y \(x_i = (x_{i1}, \ldots, x_{ip})^T \in \mathbb{R}^p\). El modelo de regresión lineal asume que
\[ Y_i = r(x_i) + \varepsilon_i \equiv \sum_{j=1}^{p} \beta_j x_{ij} + \varepsilon_i, \quad i = 1, \ldots, n, (*) \]
donde \(\mathbb{E}(\varepsilon_i) = 0\) y \(\mathbb{V}(\varepsilon_i) = \sigma^2\).
Usualmente, queremos incluir una ordenada al origen en el modelo, por lo que adoptaremos la convención de que \(x_{i1} = 1\).
La matriz de diseño \(X\) es la matriz \(n \times p\) definida por
\[ X = \begin{pmatrix} x_{11} & x_{12} & \cdots & x_{1p} \\ x_{21} & x_{22} & \cdots & x_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ x_{n1} & x_{n2} & \cdots & x_{np} \end{pmatrix}. \]
El conjunto \(\mathcal{L}\) de vectores que se pueden obtener como combinaciones lineales de las columnas de \(X\), se denomina el espacio columna de \(X\).
Sea \(Y = (Y_1, \ldots, Y_n)^T\), \(\varepsilon = (\varepsilon_1, \ldots, \varepsilon_n)^T\), y \(\beta = (\beta_1, \ldots, \beta_p)^T\). Podemos entonces escribir (*) como
\[ Y = X\beta + \varepsilon. \tag{5.5} \]
El estimador de mínimos cuadrados \(\hat{\beta} = (\hat{\beta}_1, \ldots, \hat{\beta}_p)^T\) es el vector que minimiza la suma de cuadrados de los residuos
\[ \text{RSS} = (Y - X\beta)^T (Y - X\beta) = \sum_{i=1}^{n} \left( Y_i - \sum_{j=1}^{p} x_{ij} \beta_j \right)^2. \]
Si el problema considera un predictor \(q\)-dimensional y existe \((\mathbf{X}^T\mathbf{X})^{-1}\), entonces el estimador de mínimos cuadrados \(\hat{\boldsymbol{\beta}}\) del vector de parámetros \(\boldsymbol{\beta}\) está dado por:
\[ \hat{\boldsymbol{\beta}} = (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \mathbf{y} \]
donde \(\mathbf{X}\) y \(\mathbf{y}\) han sido definidos previamente para el modelo lineal general bajo el supuesto \(\mathbb{E}(\boldsymbol{\varepsilon}) = 0\), \(\hat{\boldsymbol{\beta}}\) es un estimador insesgado de \(\boldsymbol{\beta}\).
Por lo tanto para los regresores \(x\), \(x^2\), tenemos que:
set.seed(123)
x <- seq(1, 10, length.out = 100)
y <- exp(0.3 * x) + rgamma(100, shape = 2, rate = 0.5)
# Ajuste lineal
model_cuadrado <- lm(y ~ poly(x,defree=2))
df <- data.frame(x = x, y = y)
ggplot(df, aes(x, y)) +
geom_point() +
geom_line(aes(y = predict(model_cuadrado)), color = "red") +
ggtitle("Ajuste lineal con un polinomio de grado 2")
summary(model_cuadrado)#>
#> Call:
#> lm(formula = y ~ poly(x, defree = 2))
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -3.8361 -1.7778 -0.1533 1.0525 7.6576
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 10.4204 0.2228 46.78 < 2e-16 ***
#> poly(x, defree = 2)1 52.7816 2.2275 23.70 < 2e-16 ***
#> poly(x, defree = 2)2 19.9797 2.2275 8.97 2.25e-14 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 2.228 on 97 degrees of freedom
#> Multiple R-squared: 0.8687, Adjusted R-squared: 0.866
#> F-statistic: 321 on 2 and 97 DF, p-value: < 2.2e-16Así tenemos un modelo que para ser mas adecuado, no obstante:
Lo que se confirma con:
ad_test <- ad.test(model_cuadrado$residuals)
ad_test#>
#> Anderson-Darling normality test
#>
#> data: model_cuadrado$residuals
#> A = 1.8024, p-value = 0.0001216Supuestos Clásicos de Regresión Lineal
La regresión lineal convencional se basa en dos pilares fundamentales:
Consecuencias de Violar los Supuestos
Cuando se ignoran estos supuestos, los procedimientos inferenciales tradicionales pierden validez, pues las propiedades óptimas de los estimadores, obtenidas mediante derivación de la función de verosimilitud \(L\), dependen críticamente de la normalidad
ad_test <- ad.test(model$residuals)
ad_test#>
#> Anderson-Darling normality test
#>
#> data: model$residuals
#> A = 1.3389, p-value = 0.001704| Estadístico_AD | Valor_p | |
|---|---|---|
| A | 1.338856 | 0.0017036 |
En base a la Figura@ref(fig:supuestos-1),@ref(fig:supuestos-2),@ref(fig:supuestos-3), se observa que los supuestos del modelo no se cumplen adecuadamente, especialmente la normalidad de los residuos (como lo indica el test de Anderson-Darling) (Tabla @ref(tab:prueba-resultado)) Por lo tanto, se opta por aplicar una transformación Box-Cox, con el objetivo de mejorar el ajuste del modelo y aproximarse al cumplimiento de los supuestos.
Así aplicando el \(\lambda\) optimo, tenemos lo siguiente:
bc <- boxcox(y ~ x)
Por lo tanto nuestro nuevo modelo:
summary(lm_transformed)#>
#> Call:
#> lm(formula = y_transformed ~ x)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -1.04491 -0.38752 -0.06783 0.37187 1.99376
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 0.96133 0.14298 6.724 1.19e-09 ***
#> x 0.43473 0.02346 18.529 < 2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 0.6157 on 98 degrees of freedom
#> Multiple R-squared: 0.7779, Adjusted R-squared: 0.7757
#> F-statistic: 343.3 on 1 and 98 DF, p-value: < 2.2e-16| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 0.961 | 0.143 | 6.724 | 1.19e-09 |
| x | 0.435 | 0.023 | 18.529 | < 2e-16 |
Así podemos notar que por (Tabla @ref(tab:summary-modelo-BC)) ambos coeficientes son significantes para el modelo, además:
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
p1 <- ggplot(data.frame(x, y), aes(x, y)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE, color = "red") +
ggtitle("Datos Originales")
p2 <- ggplot(data.frame(x, y_transformed), aes(x, y_transformed)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE, color = "blue") +
ggtitle(paste("Datos Transformados (λ =", round(lambda, 2), ")"))
ad_test <- ad.test(lm_transformed$residuals)
ad_test#>
#> Anderson-Darling normality test
#>
#> data: lm_transformed$residuals
#> A = 0.83135, p-value = 0.03093| Estadístico_AD | Valor_p | |
|---|---|---|
| A | 0.8313465 | 0.0309288 |
(Tabla @ref(tab:prueba-resultado))
No obstante vemos que por las Figuras (@ref(fig:supuestos-BC-1),@ref(fig:supuestos-BC-2),@ref(fig:supuestos-BC-3)) y en especial a la (Figura, @ref(fig:supuestos-BC-2)), el supuesto de normalidad no es razonable, lo cual se confirma gracias a (Tabla @ref(tab:prueba-resultado-BC))
Finalmente para el módelo multiple:
lambda <- bc$x[which.max(bc$y)]
y_transformed <- (y^lambda - 1)/lambda
# Modelo transformado
model_bc <- lm(y_transformed ~ poly(x, degree = 2))
# Gráfico del ajuste
ggplot(data.frame(x = x, y_transformed = y_transformed), aes(x, y_transformed)) +
geom_point() +
geom_line(aes(y = predict(model_bc)), color = "red", linewidth = 1) +
ggtitle("Ajuste polinómico con transformación Box-Cox ") +
theme_minimal()
Notemos lo siguiente:
La regresión lineal ha servido como base fundamental para el desarrollo de métodos no paramétricos, en el enfoque lineal tradicional:
En contraste, la regresión no paramétrica:
Métodos
El análisis de regresión no paramétrica comparte los objetivos de su contraparte paramétrica: estimar y contrastar propiedades de la función de regresión que describe la relación entre variables. Sin embargo:
\[ y = f(x) + \varepsilon \]
La única condición sobre \(\varepsilon\) es que sean v.a.i.i.d con medía 0 y varianza \(\sigma^2\)
La función de regresión \(f\) es desconocida y no se impone una forma funcional específica
Se requieren supuestos mínimos sobre \(f\):
Estas condiciones de suavidad permiten controlar el comportamiento local de las estimaciones.
Suavización
Supongamos que disponemos de una variable explicativa \(X\) y de una respuesta \(Y\) y que el modelo es aplicable al par \((X, Y )\), así sea f una función perteneciente a algún espacio de funciones \(\Omega\), con dicho espacio dependiendo de que tan suave es razonable asumir a \(f\),para estimar puntualmente a \(f\) en \(x\): \[\hat{f}(x)=\sum_{i=1}^nK(x,x_i;\lambda) y_i \] Es decir una suma ponderada donde \(K(x,x_i;\lambda)\) son pesos que dependen de \(\{x_i, i=1,..,n\}\) y del parámetro \(\lambda\) que es definido por el diseñador.
La siguiente imagen representa la probabilidad \(\mathbb{P}(x - h < X < x + h)\) y su aproximación mediante una densidad:
IMAGEN
El estimador de densidad de núcleo se basa en una aproximación local vea imagen de la densidad \(f(x)\) utilizando los datos muestrales \(X_1, \dots, X_n \sim F\) (i.i.d.).
La idea fundamental es que la probabilidad de que una observación \(X\) caiga en un intervalo pequeño alrededor de un punto \(x\) puede aproximarse, para algún \(h > 0\) conocido como parámetro de suavizamiento o ancho de banda.
Para estimar esta probabilidad usando los datos disponibles, se puede contar la proporción de observaciones que caen dentro del intervalo \([x - h, x + h]\): Cuando el ancho de banda \(h\) es pequeño, podemos aproximar la integral usando el valor de \(f\) en el punto central \(x\):
\[ \int_{x-h}^{x+h} f(t) \, dt \approx 2h \cdot f(x) \] Pues por el Teorema del Valor Medio para integrales, existe \(c \in [x-h, x+h]\) tal que: \[ \int_{x-h}^{x+h} f(t) \, dt = f(c) \cdot 2h \] Cuando \(h \to 0\), \(c \to x\).
Usando la función de distribución empírica \(F_n\), estimamos la probabilidad como la proporción de puntos muestrales en el intervalo:
\[ \mathbb{P}(x - h \leq X \leq x + h) \approx \frac{1}{n} \sum_{i=1}^n \mathbf{1}_{\{|X_i - x| \leq h\}} \]
Igualando las expresiones:
\[ 2h \cdot f(x) \approx \frac{1}{n} \sum_{i=1}^n \mathbf{1}_{\{|X_i - x| \leq h\}} \]
Despejando \(f(x)\) obtenemos el estimador inicial:
\[ f(x) \approx \frac{1}{2h} \cdot \frac{1}{n} \sum_{i=1}^n \mathbf{1}_{\{|X_i - x| \leq h\}} = \hat{f}_h(x) \] Este estimador puede verse como una versión suavizada de un histograma centrado en \(x\). Para mejorar la suavidad, se reemplaza la función indicadora \(\mathbf{1}_{\{|X_i - x| \leq h\}}\) por alguna otra una función , típicamente una función simétrica y suave como la función gaussiana.
Así se obtiene el :
\[ \hat{f}_h(x) = \frac{1}{nh} \sum_{i=1}^n K\left( \frac{x - X_i}{h} \right). \]
Este estimador distribuye la “masa” de cada observación \(X_i\) a lo largo del eje usando el núcleo \(K\), centrado en \(X_i\) y escalado por \(h\), el resultado es una función suavizada que estima la densidad \(f(x)\).
Para el parámetro \(h\)
Si los pesos \(K(x, xi; λ)\), \(i = 1, . . . , n\) son tales que:
De lo contrario para los valores **dentro de la vecindad* se incluyen constante \(c_i\) tal que:
-\(\sum_{i=1}^n c_iK(x, x_i; λ)=1\)
Así notemos lo siguiente, para estimar una densidad, necesitamos de lo siguiente:
¿Cual es la configuración optima?
set.seed(123)
datos_unif <- runif(5000, min = -3, max = -1)
### Función de densidad teórica
densidad_teorica_unif <- function(x) {
dunif(x, min = -3, max = -1)
}
### Kernel rectangular
kernel_rectangular <- function(x, centro, h) {
ifelse(abs(x - centro) <= h, 1/(2*h), 0)
}
### Estimador de densidad
estimador_densidad <- function(x, datos, h) {
mean(kernel_rectangular(x, datos, h))
}
### Rejilla de valores
rejilla <- seq(-4, 0, length.out = 500)
h1 <- 0.1
estimada_unif_h1 <- sapply(rejilla, estimador_densidad, datos = datos_unif, h = h1)
h2 <- 0.001
estimada_unif_h2 <- sapply(rejilla, estimador_densidad, datos = datos_unif, h = h2)
par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
### Gráfico 1: h = 0.2
plot(rejilla, densidad_teorica_unif(rejilla), type = "l", lwd = 3, col = "red",
main = paste("Kernel Rectangular (h =", h1, ")"),
xlab = "x", ylab = "Densidad", ylim = c(0, 0.6))
lines(rejilla, estimada_unif_h1, lwd = 2, col = "blue")
rug(datos_unif, col = "gray50")
legend("topright",
legend = c("Teórica U(-3,-1)", paste("Kernel (h =", h1, ")")),
col = c("red", "blue"), lwd = c(3, 2), cex = 0.8)
### Gráfico 2: h = 0.05
plot(rejilla, densidad_teorica_unif(rejilla), type = "l", lwd = 3, col = "red",
main = paste("Kernel Rectangular (h =", h2, ")"),
xlab = "x", ylab = "Densidad", ylim = c(0, 0.6))
lines(rejilla, estimada_unif_h2, lwd = 2, col = "blue")
rug(datos_unif, col = "gray50")
legend("topright",
legend = c("Teórica U(-3,-1)", paste("Kernel (h =", h2, ")")),
col = c("red", "blue"), lwd = c(3, 2), cex = 0.8)
set.seed(123)
datos_unif <- runif(50, min = -3, max = -1)
### Función de densidad teórica
densidad_teorica_unif <- function(x) {
dunif(x, min = -3, max = -1)
}
### Kernel rectangular
kernel_rectangular <- function(x, centro, h) {
ifelse(abs(x - centro) <= h, 1/(2*h), 0)
}
### Estimador de densidad
estimador_densidad <- function(x, datos, h) {
mean(kernel_rectangular(x, datos, h))
}
### Rejilla de valores
rejilla <- seq(-4, 0, length.out = 500)
h1 <- 0.9
estimada_unif_h1 <- sapply(rejilla, estimador_densidad, datos = datos_unif, h = h1)
h2 <- 0.05
estimada_unif_h2 <- sapply(rejilla, estimador_densidad, datos = datos_unif, h = h2)
par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
### Gráfico 1: h = 0.2
plot(rejilla, densidad_teorica_unif(rejilla), type = "l", lwd = 3, col = "red",
main = paste("Kernel Rectangular (h =", h1, ")"),
xlab = "x", ylab = "Densidad", ylim = c(0, 0.6))
lines(rejilla, estimada_unif_h1, lwd = 2, col = "blue")
rug(datos_unif, col = "gray50")
legend("topright",
legend = c("Teórica U(-3,-1)", paste("Kernel (h =", h1, ")")),
col = c("red", "blue"), lwd = c(3, 2), cex = 0.8)
### Gráfico 2: h = 0.05
plot(rejilla, densidad_teorica_unif(rejilla), type = "l", lwd = 3, col = "red",
main = paste("Kernel Rectangular (h =", h2, ")"),
xlab = "x", ylab = "Densidad", ylim = c(0, 0.6))
lines(rejilla, estimada_unif_h2, lwd = 2, col = "blue")
rug(datos_unif, col = "gray50")
legend("topright",
legend = c("Teórica U(-3,-1)", paste("Kernel (h =", h2, ")")),
col = c("red", "blue"), lwd = c(3, 2), cex = 0.8)
El riesgo cuadrático medio (\(R(\hat{f})\)) de un estimador \(\hat{f}\) se define como: \[ R(\hat{f}) = \mathbb{E}\left[(f(x) - \hat{f}(x))^2\right]. \]
Bajo condiciones de regularidad, puede descomponerse en: \[ R(\hat{f}) = \underbrace{\left(\mathbb{E}[\hat{f}(x)] - f(x)\right)^2}_{\text{Sesgo al cuadrado } (B^2)} + \underbrace{\text{Var}(\hat{f}(x))}_{\text{Varianza}}. \]
Demostración:
\[\begin{align*} R(\hat{f}) &= \mathbb{E}[(f(x) - \hat{f}(x))^2] \\ &= \mathbb{E}[(f(x) - \mathbb{E}[\hat{f}(x)] + \mathbb{E}[\hat{f}(x)] - \hat{f}(x))^2] \\ &= \underbrace{(f(x) - \mathbb{E}[\hat{f}(x)])^2}_{B^2} + \underbrace{\mathbb{E}[(\hat{f}(x) - \mathbb{E}[\hat{f}(x)])^2]}_{\text{Var}} + 2 {\mathbb{E}[(f(x) - \mathbb{E}[\hat{f}(x)])(\mathbb{E}[\hat{f}(x)] - \hat{f}(x))]}.\\ &= \underbrace{(f(x) - \mathbb{E}[\hat{f}(x)])^2}_{B^2} + \underbrace{\mathbb{E}[(\hat{f}(x) - \mathbb{E}[\hat{f}(x)])^2]}_{\text{Var}} \end{align*}\]
set.seed(123)
datos <- runif(50, -3, -1)
rejilla <- seq(-3.5, -0.5, length.out = 500)
kernel_rect <- function(x, c, h) ifelse(abs(x-c) <= h, 1/(2*h), 0)
estimador <- function(x, datos, h) rowMeans(sapply(datos, kernel_rect, x=x, h=h))
### riesgo
riesgo <- function(h, datos, grid) {
est <- estimador(grid, datos, h)
teor <- dunif(grid, -3, -1)
mse <- mean((est - teor)^2) + mean(apply(sapply(datos, kernel_rect, x=grid, h=h), 1, var))/length(datos)
return(mse)
}
### Optimización
h_opt <- optimize(riesgo, c(0.05, 1.5), datos=datos, grid=rejilla)$minimum
cat("h óptimo:", round(h_opt, 3), "\n")#> h óptimo: 0.426
Confirmando el resultado, tenemos la Regla de Silverman \[ h_{opt}\approx n^{-1/5} \]
Que en nuestra caso sería, para \(50\) datos \[ h_{opt} \approx 0.4573\]
Propiedades kernel
Si ahora le pedimos a \(K: \mathbb{R} \to \mathbb{R}^+\) sea una función continua que cumple:
De esta manera se \(f_n(x)\) serán una densidad, pues: \[\begin{align*} \int_{-\infty}^{\infty}\hat{f_h} (x)&= \int_{-\infty}^{\infty} \frac{1}{nh} \sum_{i=1}^n K\left( \frac{x - X_i}{h} \right)dx\\ &=\dfrac{1}{nh}\sum_{i=1}^n\int_{-\infty}^{\infty}K\left( \frac{x - X_i}{h} \right)dx\\ &=\dfrac{1}{nh}\sum_{i=1}^n h\int_{-\infty}^{\infty}K\left(u\right)du\\ &=\dfrac{1}{nh}\sum_{i=1}^n h=1 \\ \end{align*}\]
Así en base a Berrendero (2023) para evaluar la calidad del estimador de densidad \(\hat{f}(x)\) en un punto fijo \(x\), se consideran dos medidas fundamentales:
Sesgo:
\[ \text{Sesgo}[\hat{f}(x)] = E[\hat{f}(x)] - f(x) \]
Varianza:
\[ \text{Var}[\hat{f}(x)] = E\left[(\hat{f}(x) - E[\hat{f}(x)])^2\right] \]
La consistencia débil del estimador \(\hat{f}(x)\) requiere que:
De esta manera tenemos que :
Así notemos lo siguiente ,el estimador kernel de densidad es:
\[ \hat{f}_h(x) = \frac{1}{nh}\sum_{i=1}^n K\left(\frac{x-X_i}{h}\right) \]
Como \(X_i\) son \(v.a.i.i.d\): \[\begin{align*} \mathbb{E}[\hat{f}_h(x)] &= \frac{1}{h}\mathbb{E}\left[K\left(\frac{x-X}{h}\right)\right] \\ &= \frac{1}{h}\int K\left(\frac{x-t}{h}\right)f(t)dt \\ &\text{Con el cambio de variable } u = \frac{x-t}{h}: \\ &= \int K(u)f(x-hu)du \end{align*}\]
Usando expansión de Taylor para \(h \to 0\):
\[\begin{align*} f(x-hu) &\approx f(x) - huf'(x) + \mathcal{O}(h^2) \\ \mathbb{E}[\hat{f}_h(x)] &\approx f(x)\underbrace{\int K(u)du}_{=1} - hf'(x)\underbrace{\int uK(u)du}_{=0} + \mathcal{O}(h^2) \\ &\to f(x) \text{ cuando } h \to 0 \end{align*}\]
Por otro lado la varianza, dado que \(X_i\) son \(v.a.i.i.d\)
\[\begin{align*} \text{Var}[\hat{f}_h(x)] &= \text{Var}\left[\frac{1}{nh}\sum_{i=1}^n K\left(\frac{x-X_i}{h}\right)\right] \\ &= \frac{1}{n^2h^2}\sum_{i=1}^n \text{Var}\left[K\left(\frac{x-X_i}{h}\right)\right] \quad \\ &= \frac{1}{nh^2}\text{Var}\left[K\left(\frac{x-X}{h}\right)\right] \end{align*}\]
Desarrollamos ahora la varianza del término kernel:
\[\begin{align*} \text{Var}\left[K\left(\frac{x-X}{h}\right)\right] &= \mathbb{E}\left[K^2\left(\frac{x-X}{h}\right)\right]_A - \left(\mathbb{E}\left[K\left(\frac{x-X}{h}\right)\right]\right)^2_B \\ &= \underbrace{\int K^2\left(\frac{x-t}{h}\right)f(t)dt}_{\text{A}} - \underbrace{\left(\int K\left(\frac{x-t}{h}\right)f(t)dt\right)^2}_{\text{ B}} \end{align*}\]
Desarrollo del Término A (cambio de variable u = (x-t)/h):
\[\begin{align*} \int K^2\left(\frac{x-t}{h}\right)f(t)dt &= h\int K^2(u)f(x-hu)du \\ &\approx hf(x)\int K^2(u)du \quad \text{(para h pequeño, f(x-hu) ≈ f(x))} \end{align*}\]
Desarrollo del Término B :
\[\begin{align*} \left(\int K\left(\frac{x-t}{h}\right)f(t)dt\right)^2 &= \left(h\int K(u)f(x-hu)du\right)^2 \\ &\approx \left(hf(x)\int K(u)du\right)^2 \\ &= h^2f^2(x) \end{align*}\] \[\begin{align*} \text{Var}[\hat{f}_h(x)] &\approx \frac{1}{nh^2}\left[hf(x)\int K^2(u)du - h^2f^2(x)\right] \\ &= \frac{f(x)}{nh}\int K^2(u)du - \frac{f^2(x)}{n} \end{align*}\]
Para h → 0 y nh → ∞, el primer término domina (es de orden O(1/(nh)) mientras que el segundo es O(1/n), por lo que:
\[\begin{align} \text{Var}[\hat{f}_h(x)] &\approx \frac{f(x)}{nh}\int K^2(u)du \quad \text{(aproximación asintótica)} \end{align}\]
Estimador con series de Fourier
Para establecer un marco general, asumiremos que disponemos de observaciones de la variable de respuesta \(Y\) para \(n\) valores predeterminados de una variable independiente \(X\). Las \(n\) observaciones bivariadas disponibles, denotadas \((x_1, y_1), \ldots, (x_n, y_n)\), siguen el modelo
\[\begin{equation} \label{eq:regresion} y_i = f(x_i) + \epsilon_i, \quad i = 1, \ldots, n (\#eq:Reglin) \end{equation}\]donde \(\epsilon = (\epsilon_1, \ldots, \epsilon_n)^T\) es un vector de errores aleatorios no correlacionados que tienen media cero y varianza común \(\sigma^2\), \(f\) es una función de regresión desconocida y se satisface que \(0 \leq x_1 < \cdots < x_n \leq 1\).
Nuestro propósito es estimar \(f\), para lo cual buscaremos construir estimadores lineales que puedan escribirse en la siguiente forma general, que para un \(\lambda\) dado es una combinación lineal de las observaciones \(y_i\), donde \(K(\cdot, x_i; \lambda), \, i = 1, \ldots, n\) es una colección de funciones ponderadoras que dependen de los \(x_i\) y de un parámetro de suavización denotado\(\lambda\):
\[\begin{equation} f_\lambda(x)=\sum_{i=1}^nK(x, x_i; \lambda)y_i \end{equation}\]
Al final del siglo XVIII, Fourier propuso un método para aproximar funciones periódicas usando combinaciones lineales de funciones trigonométricas sencillas. Sin embargo, aunque el método propuesto por Fourier se refiere a funciones continuas, sus ideas se pueden aplicar mejor en el contexto más amplio de las funciones cuadrado integrables.
Definición: EL espacio \(L_2[0,1]\) se define como el espacio funcional de todas las funciones cuadrado integrables en [0, 1], es decir \[L_2[0,1]=\{f:[0,1] \rightarrow \mathbb{R}|\int_{0}^1 |f(x)|^2dx<\infty\}\]
Asumamos entonces que \(f \in L_2[0, 1]\) y que puede representarse usando la expansión: \[\begin{equation} f(x) = \sum_{j=1}^{\infty} \beta_j f_j(x) (\#eq:Fourier1) \end{equation}\]
donde los \(\beta_j\), \(j = 1, 2, \ldots\) son los , definidos en el Anexo @ref(sec-apendice), y las funciones \(f_j\), \(j = 1, 2, \ldots\) conforman una colección de funciones en \(L_2[0, 1]\). Diremos que la colección de funciones \(\{f_j\}_{j=1}^\infty\) utilizada para construir la expansión @ref(eq:Fourier1) es una (CONS), la definición de CONS se encuentra en el Anexo @ref(sec-apendice).
En consecuencia, el modelo @ref(eq:Reglin) puede representarse como:
\[\begin{equation} y_i = \sum_{j=1}^{\infty} \beta_j f_j(x_i) + \varepsilon_i, \quad i = 1, 2, \ldots, n \end{equation}\]
Esto significa que los datos siguen un modelo lineal con una cantidad infinita de coeficientes de regresión desconocidos. Ahora bien, un resultado conocido sobre los coeficientes de Fourier, llamado la Relación de Parseval, definida como \[ \sum_{j=1}^{\infty} \beta_j^2 = \|f\|^2. \] Esta relación nos permite concluir que los coeficientes \(\beta_j\) deben decaer a cero eventualmente. Entonces uno podría entonces asumir que existe un entero \(\lambda\) tal que:
\[ f \approx \sum_{j=1}^{\lambda} \beta_j f_j \] y, por tanto, podríamos escribir la siguiente aproximación:
\[\begin{equation} y_i \approx \sum_{j=1}^{\lambda} \beta_j f_j(x_i) + \varepsilon_i, \quad i = 1, 2, (\#eq:flam) \end{equation}\]
Este modelo @ref(eq:flam) tiene la forma de un modelo lineal, por lo que una posible solución al problema de la estimación de \(f\) sería estimar los coeficientes de Fourier \(\{\beta_j\}_{j=1}^{\lambda}\) utilizando la técnica de mínimos cuadrados.
Para ello, definimos la matriz:
\[ X_\lambda = \{f_j(x_i)\}_{i=1,\ldots,n;\; j=1,\ldots,\lambda} \]
Si la matriz \((X_\lambda^\top X_\lambda)^{-1}\) existe, entonces el estimador de \(\boldsymbol{\beta}\) tendrá la forma general:
\[\begin{equation} \boldsymbol{\beta}_\lambda = (X_\lambda^\top X_\lambda)^{-1} X_\lambda^\top \mathbf{y} (\#eq:minlam) \end{equation}\]
con \(\mathbf{y} = (y_1, y_2, \ldots, y_n)^\top\). Entonces, nuestro estimador de \(f\) será:
\[\begin{equation} f_\lambda(x) = \sum_{j=1}^{\lambda} \beta_{\lambda j} f_j(x) (\#eq:estlam) \end{equation}\]
Supongamos que decidimos utilizar este estimador de \(f\). Quedarían en ese caso algunas dudas por resolver, recordando que \(\sum_{j=1}^{\lambda} \beta_j f_j\) es solamente una aproximación a \(f\) y que la verdadera función de regresión es, en realidad, de dimensión infinita.
Estimacion de la varianza {#subsec-varianzas}
Una vez estimada la función de regresión, el análisis de regresión requiere estimar la varianza \(\sigma^2\) para poder continuar con el proceso de inferencia estadística. Un primer acercamiento podría ser adaptar el estimador de varianza del análisis de regresión lineal:
\[\begin{equation} \sigma^2_\lambda = \frac{\text{SCE}_\lambda}{n - \lambda} (\#eq:sigma1) \end{equation}\]
donde la suma de cuadrados de los residuales \(\text{SCE}_\lambda\) se define como:
\[\begin{equation} \text{SCE}_\lambda = \sum_{i=1}^n \left(y_i - f_\lambda(x_i)\right)^2 \end{equation}\]
Sin embargo, esta definición tiene algunas dificultades, porque el estimador \(\sigma^2_\lambda\) definido en @ref(eq:sigma1) depende del valor de \(\lambda\), que ya se ha utilizado para estimar \(f\). A diferencia del estimador de varianza en el modelo de regresión lineal, el estimador \(\sigma^2_\lambda\) no es, en general, insesgado. El sesgo del estimador \(\sigma^2_\lambda\) es:
\[\begin{equation} B(\sigma^2_\lambda) = \frac{f^\top (I - S_\lambda) f}{n - \lambda} \end{equation}\]
donde \(f = \big(f(x_1), f(x_2), \ldots, f(x_n)\big)^\top\) y la matriz de suavizamiento es:
\[ S_\lambda = X_\lambda (X_\lambda^\top X_\lambda)^{-1} X_\lambda^\top \]
Estudios previos concluyen que una elección de \(\lambda\) que minimice el sesgo del estimador \(\sigma^2_\lambda\) no necesariamente conduce a una buena selección del estimador \(f_\lambda\) de \(f\).
Una estimación de la varianza \(\sigma^2\) que no dependa de \(\lambda\) permitiría hacerla independiente del estimador de \(f\). En este contexto de los modelos no paramétricos se han propuesto varios estimadores de la varianza que no dependen de la elección de \(\lambda\). Algunos estimadores son:
El estimador HKT es un estimador insesgado de \(\sigma^2\) que generaliza los dos estimadores previos de Rice y GSJS. El estimador de Rice coincide con el estimador HKT cuando \(r = 1\), con \(m_1 = 0\) y \(m_2 = 1\). Y el estimador GSJS se obtiene a partir del estimador HKT con \(m_1 = 0\) y \(m_2 = 2\), si el diseño es equidistante.
Algunas ideas sobre cómo conducir inferencias
Tal como en los modelos lineales, no parece razonable conservar una función \(f_j\) que tenga asociado un \(\beta_j = 0\). Así que una primera prueba formal del análisis de regresión será contrastar las hipótesis \[ H_0 : \beta_j = 0\quad\text{vs}\quad H_1 : \beta_j \neq 0 \] Un posible estadístico de prueba, adaptado de los modelos lineales, sería: \[ Z_j = \frac{\beta_j}{\sqrt{(X^\top_\lambda X_\lambda)_j \sigma^2}} \] donde \((X^\top_\lambda X_\lambda)_j\) es el elemento \(j\) de la diagonal principal de la matriz \(X^\top_\lambda X_\lambda\). En lugar de \(\sigma^2\) podríamos utilizar uno de los estimadores propuestos en la Sección @ref(subsec-varianzas) De igual forma, un intervalo de confianza para \(f\) sería: \[ f_\lambda(x) \pm Z_{\alpha/2} \sqrt{\sigma^2 f_\lambda(x)^\top (X^\top_\lambda X_\lambda)^{-1} f_\lambda(x)} \] con \(f_\lambda(x) = (f_1(x), \dots, f_\lambda(x))^\top\). El análisis de homocedasticidad y demás pruebas de diagnóstico se podrían tomar prestados de los métodos clásicos de análisis de regresión.
Ejemplo El siguiente ejemplo fue tomado de Univalle, los datos representados en la siguinete tabla corresponden a mediciones del promedio horario de Monóxido de Carbono (CO) en una estación ubicada en el centro de la ciudad de Cali, Colombia (Calle 15), el 1 de marzo de 2004.
| HORA | CO | HORA | CO |
|---|---|---|---|
| 0 | 2.24 | 12 | 3.22 |
| 1 | 1.53 | 13 | 2.94 |
| 2 | 1.80 | 14 | 2.13 |
| 3 | 0.87 | 15 | 2.88 |
| 4 | 1.21 | 16 | 1.95 |
| 5 | 1.02 | 17 | 1.51 |
| 6 | 1.39 | 18 | 1.54 |
| 7 | 1.90 | 19 | 1.60 |
| 8 | 2.49 | 20 | 1.22 |
| 9 | 3.33 | 21 | 1.37 |
| 10 | 2.82 | 22 | 1.96 |
| 11 | 3.34 | 23 | 1.66 |
Nuestro propósito es ajustar un modelo del tipo @ref(eq:Reglin) a este conjunto de observaciones de un día, para ello hagamos un reescalamiento de las horas para estar en el intervalo [0,1].
Supongamos que disponemos de la siguiente sucesión ortonormal completa (CONS) de funciones: \[\begin{align*} f_1(x) &= 1, \\ f_j(x) &= \sqrt{2} \cos((j - 1)\pi x), \quad j = 2, 3, \dots \end{align*}\]
Usando esto, nuestro estimador @ref(eq:estlam) tomaría la siguiente forma: \[\begin{equation} f_\lambda(x) = \sum_{j=1}^{\lambda} \beta_j^\lambda \sqrt{2} \cos((j - 1)\pi x) \end{equation}\]
Podemos estimar los \(\beta_{\lambda j}\) utilizando el método de mínimos cuadrados tal como se describió en @ref(eq:minlam), supongamos que para los datos decidimos estimar \(f\) usando \(\lambda=6\), es decir queremos aproximar \(f\) con las primeras 6 funciones de la CONS de cosenos:
\[\begin{equation} f_6(x) = \beta_{61} f_1 + \beta_{62} f_2 + \beta_{63} f_3 + \beta_{64} f_4+ \beta_{65} f_5+ \beta_{66} f_6 \end{equation}\]
# Datos
hora <- 0:23
co <- c(2.24, 1.53, 1.80, 0.87, 1.21, 1.02, 1.39, 1.90,
2.49, 3.33, 2.82, 3.34, 3.22, 2.94, 2.13, 2.88,
1.95, 1.51, 1.54, 1.60, 1.22, 1.37, 1.96, 1.66)
# Reescalar las horas al intervalo [0,1]
hora_reescalada <- (hora - min(hora)) / (max(hora) - min(hora))
# Definir las funciones f_j(x) para j = 2,...,6
f2 <- function(x) sqrt(2) * cos(1 * pi * x)
f3 <- function(x) sqrt(2) * cos(2 * pi * x)
f4 <- function(x) sqrt(2) * cos(3 * pi * x)
f5 <- function(x) sqrt(2) * cos(4 * pi * x)
f6 <- function(x) sqrt(2) * cos(5 * pi * x)
# Matriz X_\lambda
matriz <- cbind(
1, # Columna de unos (intercepto)
f2(hora_reescalada),
f3(hora_reescalada),
f4(hora_reescalada),
f5(hora_reescalada),
f6(hora_reescalada)
)
# Beta
beta <- solve(t(matriz) %*% matriz) %*% t(matriz) %*% co
# Mostrar los estimadores
print(beta)#> [,1]
#> [1,] 2.007396347
#> [2,] -0.032968097
#> [3,] -0.562119784
#> [4,] 0.002013006
#> [5,] 0.380031077
#> [6,] 0.131421828# Valores
x_nuevo <- seq(0, 1, length.out = 100)
# Estimación de f
funcion_suma <- function(x) {
beta[1] +
beta[2] * sqrt(2) * cos(1 * pi * x) + # f2
beta[3] * sqrt(2) * cos(2 * pi * x) + # f3
beta[4] * sqrt(2) * cos(3 * pi * x) + # f4
beta[5] * sqrt(2) * cos(4 * pi * x) + # f5
beta[6] * sqrt(2) * cos(5 * pi * x) # f6
}
# Evaluar en los nuevos puntos
y_ajustado <- sapply(x_nuevo, funcion_suma)
# Graficar datos con la función ajustada
plot(hora_reescalada, co,
main = "Ajuste de funciones base coseno",
xlab = "Hora (reescalada a [0,1])",
ylab = "Concentración de CO",
pch = 19, col = "blue", cex = 0.7)
lines(x_nuevo, y_ajustado, col = "red", lwd = 2)
El modelo de regresión polinómica {#subsec-polinomica}
Asumamos por un momento que disponemos de una función de regresión \(f\) conocida. De acuerdo con el Teorema de Taylor, si tenemos una función continua f en el intervalo cerrado \([a, b]\) que tiene \(\lambda+1\) derivadas continuas en el intervalo (a, b), entonces si \(x\) y \(c\) son dos puntos en (a, b), se sigue que \(f(x)\) se puede representar como: \[\begin{equation} f(x)=f(c)+(x-c)f'(c)+\frac{(x-c)^2}{2!}f''(c)+\dots+ \frac{(x-c)^\lambda}{\lambda!}f''(c)+r_\lambda(\xi). (\#eq:Taylor) \end{equation}\]
Donde el término \(r_\lambda(\xi)\) depende de la derivada \(f^{\lambda+1}\) de \(f\) para cierto \(\xi\) entre \(x\) y \(c\), y se conoce como el \(\lambda\)-residuo. La expresión de \(f\) usando el Teorema de Taylor @ref{eq:Taylor} puede reescribirse como:
\[\begin{equation} f(x)=\sum_{j=0}^\lambda\theta_jf_j(x)+r_\lambda(x). (\#eq:Taylor2) \end{equation}\]
donde \[ \theta_1 = f(c), \quad \theta_j = \frac{f^{(j-1)}(c)}{(j-1)!}, \quad j = 2, \ldots, \lambda, \quad \text{y}\quad f_j(x) = (x - c)^{j-1}, \quad j = 1, \ldots, \lambda. \]
Es decir, la función \(f\) se puede representar como la suma de un polinomio \(\pi_\lambda(x)\) de grado \(\lambda\) más un residuo \(r_\lambda(x)\).
Ahora bien, en la medida que \(r_\lambda(x)\) sea más pequeño, entonces \(\pi_\lambda(x)\) se parecerá más a \(f(x)\). En resumen, existe un \(\lambda\) tal que: \[ f(x) \approx \sum_{j=0}^{\lambda} \theta_j f_j(x) \]
Supongamos ahora que hemos colectado información sobre dos variables \(X\) y \(Y\) que se relacionan según el modelo @ref(eq:Reglin), por lo que la función de regresión no es conocida. Si asumimos que \(f\) es una función de regresión que tiene \(\lambda + 1\) derivadas continuas y puede representarse en la forma @ref(eq:Taylor2), entonces, si \(r_\lambda(x)\) es uniformemente pequeño, podríamos escribir:
\[\begin{equation} f(x) \approx \sum_{j=0}^{\lambda} \theta_j f_j(x) + \varepsilon_j (\#eq:polinomico) \end{equation}\]
Diremos entonces que nuestras observaciones se comportan aproximadamente según un modelo de regresión (no-paramétrica) polinómica, el cual presupone que los residuos \(r_\lambda(x_1), \ldots, r_\lambda(x_n)\) de la expansión de \(f\) en series de Taylor han sido incorporados en los errores aleatorios del modelo @ref(eq:polinomico).
Para encontrar los coeficientes \(\theta_j\), se podría utilizar el método de mínimos cuadrados, o cualquier otro.
Por otra parte, si \(f\) se puede aproximar bien utilizando el polinomio \(\pi_\lambda(x)\), y si los residuos en los \(x_i\) son pequeños en comparación con los errores \(\varepsilon_i\), entonces el modelo @ref(eq:polinomico) debería funcionar bien. En otro caso, podríamos estar en problemas con este modelo, porque como tanto los residuos como los errores son desconocidos.
Introducción
Para establecer un marco general, se asume que se dispone de observaciones de la variable de respuesta \(Y\) para \(n\) valores predeterminados de una variable independiente \(X\). Las \(n\) observaciones bivariadas disponibles, denotadas \((x_1, y_1), \dots, (x_n, y_n)\), siguen el modelo
\[ y_i = f(x_i) + \varepsilon_i, \quad i = 1, \dots, n (\#eq:modelo) \] donde \(\varepsilon = (\varepsilon_1, \dots, \varepsilon_n)^T\) es un vector de errores aleatorios no correlacionados que tienen media cero y varianza común \(\sigma^2\), \(f\) es una función de regresión desconocida y se satisface que \(0 \le x_1 < \cdots < x_n \le 1\).
Según Univalle, para efectos de presentación de los resultados, se asume que los valores de \(X\) se han elegidos así: \[ x_i = \frac{2i - 1}{2n}, \quad i = 1, 2, \dots, n (\#eq:particion) \] El propósito es estimar \(f\) en @ref(eq:modelo), para lo cual se busca construir estimadores lineales que puedan escribirse en la siguiente forma general, que para un \(\lambda\) dado es una combinación lineal de las observaciones \(y_i\), donde \(K(\cdot; x_i; \lambda)\), \(i = 1, \ldots, n\) es una colección de funciones ponderadoras que dependen de los \(x_i\) y de un parámetro de suavización denotado \(\lambda:\)
\[ f_\lambda(x) = \sum_{i=1}^n K(x, x_i; \lambda) y_i (\#eq:fSmooth) \]
Estimadores Kernel
De acuerdo con Univalle, se discute sobre los estimadores de series, donde se encontra que el estimador de cosenos puede representarse como un estimador lineal, que se construye con una colección de funciones de pesos, una para cada \(x_i\), que depende de \(\lambda:\)
\[ f_\lambda(x) = \frac{1}{n} \sum_{i=1}^{n} y_i K_\lambda(x; x_i) (\#eq:fCONS) \] donde:
\[ K_\lambda(x; s) = \left[1 + \cos\left(\frac{\lambda \pi (x - s)}{2}\right)\right] D\left(\frac{x - s}{2}; \lambda - 1\right) + \left[1 + \cos\left(\frac{\lambda \pi (x + s)}{2}\right)\right] D\left(\frac{x + s}{2}; \lambda - 1\right) (\#eq:kCONS) \] En @ref(eq:kCONS) la función \(D\) se conoce como kernel de Dirichlet, que proviene del análisis de series de Fourier y se define como:
La función de asignación de pesos @ref(eq:kCONS) que se utiliza en ese caso tiene la limitación de oscilar a medida que se aleja del punto de estimación \(x\), de tal manera que asigna pesos bajos a algunos puntos cercanos al punto de estimación y pesos más altos a puntos más alejados de este punto de estimación.
\[ D(u; k) = \frac{\sin(k \pi u)}{\sin(\pi u)} \]
Con este enfoque, la observación \(y_i\) asociada con un \(x_i\) cercano a \(x\) tendrá asignado un peso más alto, contribuyendo lo máximo posible a la estimación de la función de regresión en \(x\). Al mismo tiempo, observaciones \(y_i\) asociadas con \(x_i\) alejadas de \(x\) tendrán pesos cada vez más pequeños a medida que estén más alejadas de \(x\).
Una opción de construcción de estimadores que eviten esta limitación sería utilizar la misma idea pero con una función ponderadora que asigne consistentemente pesos más bajos a observaciones más alejadas del punto de estimación, al mismo tiempo que asigna pesos altos a las observaciones cercanas al punto deestimación. Una manera de producir estos estimadores es utilizar colecciones de funciones ponderadoras que tengan la forma general:
\[ K(x, x_i; \lambda) = \frac{1}{\lambda} K\left(\frac{x - x_i}{\lambda}\right), \quad i = 1, \dots, n (\#eq:KmodYEsc) \] donde \(K\) es una función simétrica con soporte en \([-1, 1]\) y su máximo en cero. Al parámetro \(\lambda\), que en los estimadores de series era el número de términos de la serie que se conservaban en la estimación de \(f\), lo llamaremos indistintamente parámetro de suavización o ancho de banda y cumplirá la misma función: servir de base para elegir el mejor estimador de \(f.\) Pero a diferencia de los estimadores de series, en los estimadores kernel \(\lambda\) no tiene que ser un entero, aunque sí se restringirá su valor a que sea cualquier real no negativo.
Con este modo de hacer, la observación \(y_i\) asociada con un \(x_i\) cercano a \(x\) tendrá asignado un peso más alto, contribuyendo de esta manera lo máximo posible a la estimación de la función de regresión en \(x\). Al mismo tiempo, observaciones \(y_i\) asociadas con \(x_i\) alejadas de \(x\) tendrán pesos cada vez más pequeños a medida que estén más alejados de \(x\). Para garantizar lo anterior, le exigiremos por lo pronto a la función \(K\) que tenga las siguientes dos propiedades:
\[ \int_{-1}^1 K(u) \, du = 1 \\ \int_{-1}^1 u K(u) \, du = 0 (\#eq:CondicK) \]
La primera propiedad en @ref(eq:CondicK) ayuda a garantizar que la suma de los pesos sea igual a 1; la segunda, a que \(K\) sea simétrica alrededor de cero.
A estas funciones se les llama funciones kernel. Y a los estimadores basados en estas funciones se les llama estimadores kernel.
Dos funciones kernel de uso común, que se conocen como:
Kernel cuadrático (Epanechnikov): \[ K(u) = \frac{3}{4}(1 - u^2) \quad \text{si } |u| \le 1, \text{ y } 0 \text{ en otro caso} \]
Kernel biponderado (biweight): \[ K(u) = \frac{15}{16}(1 - u^2)^2 \quad \text{si } |u| \le 1, \text{ y } 0 \text{ en otro caso} \]
# Definimos las funciones kernel
epanechnikov <- function(u) ifelse(abs(u) <= 1, 0.75 * (1 - u^2), 0)
biweight <- function(u) ifelse(abs(u) <= 1, (15/16) * (1 - u^2)^2, 0)
#gaussian <- function(u) dnorm(u) # sin recorte
# Valores para graficar
u_vals <- seq(-1.5, 1.5, length.out = 500)
# Graficamos
plot(u_vals, epanechnikov(u_vals), type = "l", col = "blue", lwd = 2,
ylim = c(0, 1), ylab = "K(u)", xlab = "u", main = "Funciones Kernel")
lines(u_vals, biweight(u_vals), col = "red", lwd = 2)
#lines(u_vals, gaussian(u_vals), col = "green", lwd = 2, lty = 2)
legend("topright", legend = c("Epanechnikov", "Biweight"),
col = c("blue", "red"), lty = c(1, 1, 2), lwd = 2)
Estimador kernel de Priestley-Chao
Con el estimador de cosenos se concluyo que una forma particular del estimador es:
\[ f_{\lambda}(x) = \frac{1}{n} \sum_{i=1}^n K(x; x_i; \lambda)y_i \]
Si se incorpora ahora el kernel genérico, se obtiene el estimador kernel:
\[ f_{\lambda}(x) = \frac{1}{n \lambda} \sum_{i=1}^n K\left(\frac{x - x_i}{\lambda}\right)y_i \]
Este estimador kernel es un estimador kernel genérico que luce tal como el estimador de cosenos @ref(eq:kCONS). Aunque ambos asignan pesos altos a puntos cercanos a \(x\), el kernel genérico reduce el peso más rápidamente con la distancia.
Este estimador fue propuesto por Priestley (1972) para el caso de diseños igualmente espaciados como el del diseño. Sin embargo, su desempeño no es muy bueno, como se ve en la Figura @ref(fig:KCuadraCO), donde se destaca la escasa habilidad del estimador en los extremos del diseño.
# Datos de hora y CO
hora <- 0:23
co <- c(2.24, 1.53, 1.80, 0.87, 1.21, 1.02, 1.39, 1.90, 2.49, 3.33, 2.82, 3.34,
3.22, 2.94, 2.13, 2.88, 1.95, 1.51, 1.54, 1.60, 1.22, 1.37, 1.96, 1.66)
# Normalizar la hora al intervalo [0,1]
hora_cod <- hora / max(hora)
# Kernel cuadrático (biweight)
k_biweight <- function(u) {
ifelse(abs(u) <= 1, (15/16) * (1 - u^2)^2, 0)
}
# Estimación kernel para un punto x dado
f_lambda <- function(x, x_i, y_i, lambda) {
n <- length(x_i)
sum(k_biweight((x - x_i)/lambda) * y_i) / (n * lambda)
}
# Evaluar la estimación en una malla de puntos
lambda <- 0.15
x_grid <- seq(0, 1, length.out = 200)
f_est <- sapply(x_grid, function(x) f_lambda(x, hora_cod, co, lambda))
# Graficar
plot(hora_cod, co,main="Estimación Kernel de Concentraciones de CO en Cali" ,pch = 19,col="blue" ,xlab = "HORA codificada", ylab = "CO")
lines(x_grid, f_est, type = "l")
El ancho de banda se elige de manera similar, basado en la minimización del criterio de validación cruzada generalizada GCV.
Si uno dispusiera de un segundo conjunto de \(n\) datos, estimaría \(f_\lambda\) con los primeros \(n\) datos y buscaría el \(\lambda\) que minimice \(P(\lambda)\) usando el segundo conjunto de datos. En la práctica, sin embargo, sería mejor estimar \(f\) utilizando los \(2n\) datos, lo que nos lleva a continuar nuestra búsqueda de un nuevo procedimiento para elegir \(\lambda\) que no dependa de \(\sigma^2\).
Una opción es dividir el conjunto de \(n\) observaciones en \(n\) sub-muestras de tamaño \(n - 1\) mediante el mecanismo de dejar por fuera una observación diferente cada vez. Si denotamos \(\hat{f}_\lambda^{(i)}\) a la estimación de \(f_i\) obtenida al suprimir de la muestra la observación \(i\), entonces la observación \(y_i\) sería una observación adicional que podríamos utilizar para construir un estimador de \(P(\lambda)\) que denotaremos \(\text{CV}(\lambda)\) y que llamaremos criterio de validación cruzada:
\[ \text{CV}(\lambda) = \frac{1}{n} \sum_{i=1}^n (y_i - \hat{f}_\lambda^{(i)})^2 \]
Calcular el criterio CV es complejo, porque requiere una gran carga de procesamiento computacional. Green (1993) sugieren simplificar este cálculo utilizando en cambio:
\[ \text{CV}_{GS}(\lambda) = \frac{1}{n} \sum_{i=1}^n \left( \frac{y_i - \hat{f}_{\lambda, i}}{1 - S_{\lambda, ii}} \right)^2 \]
donde \(S_{\lambda, ii}\) es el elemento \(i\) de la matriz \(S_\lambda\).
Otra posible solución se debe a Wahba (1990), quien sugiere
utilizar otro criterio llamado validación cruzada generalizada
\(\text{GCV}(\lambda)\), definido por como:
\[ \text{GCV}_{GS}(\lambda) = \frac{n^{-1} \sum_{i=1}^n (y_i - \hat{f}_{\lambda, i})^2}{(1 - n^{-1} \text{tr}(S_\lambda))^2} \]
Comparando las expresiones se ve que los residuales obtenidos al eliminar una de las observaciones se pueden obtener a partir de los residuales mediante el mecanismo de dividirlos entre el factor \(1 - S_{\lambda, ii}\). La idea del criterio GCV es reemplazar estos factores por su promedio \(1 - n^{-1} \text{tr}(S_\lambda)\).
Como en el numerador aparece la suma de cuadrados de los residuales \(\text{RSS}\), una expresión más sencilla sería Eubank (1999):
\[ \text{GCV}(\lambda) = \frac{n^{-1} \text{RSS}(\lambda)}{(n^{-1} \text{tr}(I - S_\lambda))^2} (\#eq:GCV) \]
Aunque la palabra “generalizada” deja la impresión de que el segundo criterio generaliza el primero, esto no es en general cierto y se trata de criterios diferentes que permiten estimar el riesgo de predicción \(P(\lambda)\). Wahba (1990) justifica el uso de este criterio como un buen método de selección de \(\lambda\) demostrando que si \(n^{-1} \text{tr}(S_\lambda) < 1\), entonces la diferencia entre \(\mathbb{E}[\text{GCV}(\lambda)]\) y \(P(\lambda)\), relativa al tamaño de \(R(\lambda)\), será pequeña, especialmente para tamaños de muestra grande.
En este caso el estimador lucirá como \(\hat{f}_\lambda = S_\lambda y\) , donde un término \(S_{ij}\) de la matriz \(S_\lambda\) lucirá como:
\[ S_{ij} = \frac{1}{n \lambda} K\left( \frac{x_j - x_i}{\lambda} \right), \quad i, j = 1, \dots, n \]
Otro aspecto importante es la selección del kernel. Una alternativa sería identificar el kernel que minimiza el riesgo de \(\hat{f}_h\), una vez elegido el ancho de banda. Estudios previos (Benedetti 1975, Gasser & Müller 1979) identifican el kernel cuadrático de Epanechnikov (1969), adaptado para que sea positivo en \([-1, 1]\), como una buena elección.
Otros kernels propuestos incluyen:
\[ K(u) = \frac{1}{\sqrt{2\pi}} e^{-u^2}, \quad u \in (-\infty, \infty) \]
\[ K(u) = e^{-u^2}, \quad u \in (-\infty, \infty) \]
\[ K(u) = \frac{35}{32}(1 - u^2)^3, \quad u \in [-1, 1] \]
Sin embargo, la elección del kernel tiene baja influencia comparada con la elección del ancho de banda.
Varios intentos posteriores tratan de mejorar la habilidad del estimador genérico. El primero fue propuesto por Benedetti(1975), y también por nadaraya (1964) y Watson (1964), para el caso en el que \(X\) es una variable aleatoria. Este estimador de conoce como estimador de Nadaraya-Watson y se define así:
\[ f_{\lambda}(x) = \frac{\sum_{i=1}^n K\left(\lambda^{-1}(x - x_i)\right)y_i}{\sum_{j=1}^n K\left(\lambda^{-1}(x - x_j)\right)} \]
Este estimador ajusta los pesos del kernel para que sumen 1, convirtiéndolo en un promedio ponderado.
# Normalizar la hora a [0, 1]
hora_cod <- hora / max(hora)
# Estimación kernel usando ksmooth (Nadaraya-Watson con kernel normal)
ks <- ksmooth(x = hora_cod, y = co, kernel = "normal", bandwidth = 0.15)
# Graficar resultados
plot(hora_cod, co, pch = 1, xlab = "HORA codificada", ylab = "CO")
lines(ks, col = "blue", lwd = 2)
La versión original parte de:
\[ f(x) = E(Y \mid X = x) \]
Requiriendo la distribución condicional \(h(y|x)\), y usando:
\[ f(x) = \int y h(y|x) \, dy = \int y \frac{h(x, y)}{h_X(x)} \, dy \tag{3.9} \]
Una solución es estimar densidades vía suavización. Härdle (1990b) y Simonoff (1996) discuten extensamente esta metodología.
Un estimador kernel-producto de \(h(x, y)\) es:
\[ \hat{h}(x, y) = \frac{1}{h_x h_y} \sum_{i=1}^n K_x\left(\frac{x - x_i}{h_x}\right) K_y\left(\frac{y - y_i}{h_y}\right) \]
Y un estimador de \(h_X(x)\) es:
\[\hat{h}(x) = \frac{1}{h_x} \sum_{i=1}^n K_x\left(\frac{x - x_i}{h_x}\right) \] Con estas ecuaciones se obtiene el estimador de Nadaraya-Watson.
Otro estimador kernel fue propuesto por Gasser & Müller (1979):
\[ \hat{f}(x) = \sum_{i=1}^n \left[ h^{-1} \int_{s_{i-1}}^{s_i} K\left(h^{-1}(x - s)\right) ds \right] y_i \]
con:
\[ s_0 = 0, \quad s_{i-1} \leq x_i \leq s_i, \quad i = 1, \dots, n - 1, \quad s_n = 1 \]
# Normalizar hora
x <- hora / max(hora)
y <- co
n <- length(x)
lambda <- 0.15
# Kernel normal
K <- function(u) dnorm(u)
# Crear los límites s_i
s <- numeric(n + 1)
s[1] <- x[1] - (x[2] - x[1]) / 2
s[n + 1] <- x[n] + (x[n] - x[n - 1]) / 2
for (i in 2:n) {
s[i] <- (x[i - 1] + x[i]) / 2
}
# Estimador de Gasser-Müller
gasser_muller <- function(x_eval) {
sum(sapply(1:n, function(i) {
integral <- integrate(function(s_var) K((x_eval - s_var)/lambda), lower = s[i], upper = s[i + 1])$value
(1 / lambda) * integral * y[i]
}))
}
# Evaluar en una malla de puntos
x_grid <- seq(0, 1, length.out = 200)
f_gm <- sapply(x_grid, gasser_muller)
# Graficar
plot(x, y, pch = 1, xlab = "HORA codificada", ylab = "CO")
lines(x_grid, f_gm, type = "l", col = "red", lwd = 2)
Tenemos entonces tres estimadores kernel:
Estos estimadores tienen propiedades asintóticas similares. asintóticas, por lo que en lo sucesivo no nos detendremos a identificar cuál de ellos está en uso. desigualmente espaciados, mientras que el estimador de Gasser-Muller está diseñado específicamente para este tipo de diseños desigualmente espaciados. El estimador de Nadaraya-Watson fue pensado por sus autores para diseños aleatorios y ha tomado la delantera entre la mayoría de los usuarios de los métodos de regresión no paramétrica por su flexibilidad y eficiencia.
Estimadores lineales localmente
Otro estimador común, que también es un estimador kernel, es el estimador de regresión local Cleveland (1979). Se sabe por Eubank (1999) que el estimador de Nadaraya-Watson coincide con la solución para \(\beta_0\) del problema de mínimos cuadrados ponderados:
\[ \sum_{i=1}^n (y_i - \beta_0)^2 K\left(\frac{x - x_i}{h}\right), (\#eq:RegrLoc2) \] lo que sugiere el ajuste de polinomios locales de orden \(p\), por ejemplo uno que minimice:
\[ \sum_{i=1}^n \left(y_i - \beta_0 - \dots - \beta_p(x - x_i)^p\right)^2 K\left(\frac{x - x_i}{h}\right). (\#eq:RegrLocp) \] Al estimador resultante de minimizar la suma @ref(eq:RegrLocp) se le conoce como estimador polinomial de regresión local de orden \(p\).
Implementación del estimador LOESS
Ilustraremos el caso de la solución del sistema @ref(eq:RegrLoc2), popularizado como estimador LOESS, que utiliza los \(k\) vecinos más cercanos:
La estimación LOESS de \(f(x)\) toma el valor del intercepto (\(\hat{\beta}_0\)) de la recta ajustada \(\hat{\beta}_0 + \hat{\beta}_1 x\).
# Datos
hora <- 0:23
co <- c(2.24, 1.53, 1.80, 0.87, 1.21, 1.02, 1.39, 1.90, 2.49, 3.33, 2.82, 3.34,
3.22, 2.94, 2.13, 2.88, 1.95, 1.51, 1.54, 1.60, 1.22, 1.37, 1.96, 1.66)
# Normalizar la hora a [0,1]
hora_cod <- hora / max(hora)
# Estimación LOESS (estimación lineal local con kernel Gaussiano)
modelo_loess <- loess(co ~ hora_cod, span = 0.5, degree = 2, family = "gaussian")
# Crear una malla para evaluar el modelo
x_grid <- seq(0, 1, length.out = 200)
f_loess <- predict(modelo_loess, newdata = data.frame(hora_cod = x_grid))
# Graficar
plot(hora_cod, co, pch = 1, xlab = "HORA codificada", ylab = "CO")
lines(x_grid, f_loess, col = "darkgreen", lwd = 2)
# Datos
hora <- 0:23
co <- c(2.24, 1.53, 1.80, 0.87, 1.21, 1.02, 1.39, 1.90, 2.49, 3.33, 2.82, 3.34,
3.22, 2.94, 2.13, 2.88, 1.95, 1.51, 1.54, 1.60, 1.22, 1.37, 1.96, 1.66)
# Normalizar la hora a [0,1]
hora_cod <- hora / max(hora)
# Estimación LOESS (estimación lineal local con kernel Gaussiano)
modelo_loess <- loess(co ~ hora_cod, span = 0.2, degree = 2, family = "gaussian")
# Crear una malla para evaluar el modelo
x_grid <- seq(0, 1, length.out = 200)
f_loess <- predict(modelo_loess, newdata = data.frame(hora_cod = x_grid))
# Graficar
plot(hora_cod, co, pch = 1, xlab = "HORA codificada", ylab = "CO")
lines(x_grid, f_loess, col = "darkgreen", lwd = 2)
# Datos
hora <- 0:23
co <- c(2.24, 1.53, 1.80, 0.87, 1.21, 1.02, 1.39, 1.90, 2.49, 3.33, 2.82, 3.34,
3.22, 2.94, 2.13, 2.88, 1.95, 1.51, 1.54, 1.60, 1.22, 1.37, 1.96, 1.66)
# Normalizar la hora a [0,1]
hora_cod <- hora / max(hora)
# Estimación LOESS (estimación lineal local con kernel Gaussiano)
modelo_loess <- loess(co ~ hora_cod, span = 0.5, degree = 1, family = "gaussian")
# Crear una malla para evaluar el modelo
x_grid <- seq(0, 1, length.out = 200)
f_loess <- predict(modelo_loess, newdata = data.frame(hora_cod = x_grid))
# Graficar
plot(hora_cod, co, pch = 1, xlab = "HORA codificada", ylab = "CO")
lines(x_grid, f_loess, col = "darkgreen", lwd = 2)
Estimación kernel multivariante
Todos los estimadores kernel discutidos en esta sección pueden extenderse a la estimación de funciones \(p\)-variadas. Supongamos que \(\mathbf{x} = (x_1, x_2, \dots, x_p)^T\). La solución es definir kernels multivariados como productos de kernels univariados:
\[ K(u_1, \dots, u_p) = \prod_{j=1}^p K_j(u_j), (\#eq:prodKernels) \]
El estimador de Priestley-Chao multivariado toma la forma:
\[ \hat{f}(\mathbf{x}) = \frac{1}{n} \sum_{i=1}^n y_i \prod_{h=1}^p h_h^{-1} K_h\left(h_h^{-1}(x_h - x_{ih})\right) \quad (3.16) \]
Forma general del estimador kernel multivariante
\[ \hat{f}_L(\mathbf{x}) = \frac{1}{n |L|} \sum_{i=1}^n K_p\left[L^{-1}(\mathbf{x} - \mathbf{x}_i)\right] \]
donde: - \(K_p: \mathbb{R}^p \to \mathbb{R}\) es un kernel multivariado (usualmente producto de kernels univariados). - \(L\) es una matriz \(p \times p\) de anchos de banda. - \(|L|\) es el determinante de \(L\). no-singular; y \(|L|\) es el valor absoluto del determinante de la matriz \(L\). La técnica más común para obtener \(K_p\) a partir de una función kernel univariada \(K\) es utilizar un producto de la forma @ref(eq:prodKernels).
Un ejemplo de estimación no paramétrica bivariada se presenta a continuación, en la que se estima el comportamiento del contaminante Ozono (O₃) dependiendo de la temperatura y la humedad relativa.
En este problema de la contaminación por ozono troposférico en Londres enero del 2002 , hemos ajustado un modelo bivariado de la forma:
\[ Y = f(X_1, X_2) + \varepsilon \]
donde \(Y\) representa el promedio diario de Ozono y \(X_1\) y \(X_2\) los promedios diarios de temperatura y humedad, respectivamente.
IMAGEN 3dplot.JPG
Inferencia en la estimación kernel 1. Primera condición
Garantiza que el primer término de \(\mathbb{E}[\hat{f}_h(x)]\) en (3.22) es exactamente \(f(x)\), asegurando insesgabilidad asintótica en el interior del diseño.
Segunda condición
Suprime el término de orden \(h\) en el sesgo (ver Ejercicio 2 al final del capítulo), reduciéndolo a \(\mathcal{O}(h^2 + n^{-1})\).
Demostración:
\[\int_{-1}^1 u K(u) du = 0 \implies M_1(x) = 0\]
Primera condición
Segunda condición
\(V < \infty\) garantiza varianza finita:
\[\text{var}(\hat{f}_h) \approx \frac{\sigma^2 V}{n h}\]
Estimación por intervalos
Los estimadores kernel pueden verse como promedios ponderados, por lo que el Teorema del Límite Central es aplicable y la distribución de \(f_\lambda\) será asintóticamente normal. Un posible enunciado para este resultado sería:
Proposición. Si los \(\epsilon_i\) son independientes e idénticamente distribuidos con \(E(\epsilon_i) = 0\) y \(\text{var}(\epsilon_i) = \sigma^2 < \infty\) y si \(n \to \infty\) y \(\lambda \to 0\) de tal modo que \(n\lambda \to \infty\), entonces para \(x \in [0, 1]\) se satisface que
\[ \frac{f_\lambda(x) - E[f_\lambda(x)]}{\sqrt{\text{var}[f_\lambda(x)]}} \xrightarrow{d} N[0, 1] \]
donde \(N[0, 1]\) representa una variable aleatoria normal estándar.
En el enunciado se acude únicamente al hecho de que \(f_\lambda\) es un promedio y el resultado no depende de suavización alguna, porque \(f_\lambda(x)-E[f_\lambda(x)]\) no depende de \(f\). El problema es más complejo cuando nos enfrentamos al estudio de la normalidad asintótica de la cantidad \(f_\lambda(x) - f(x)\), que sí requiere del uso de suavización. En este caso acudiremos al resultado:
\[ f_\lambda(x) - f(x) = \{f_\lambda(x) - E[f_\lambda(x)]\} + \{E[f_\lambda(x)] - f(x)\} \quad (3.25) \]
Si dividimos ambos lados de la expresión (3.25) entre \(\text{var}[f_\lambda(x)]\), entonces el problema de convergencia se reduce al segundo término de la derecha, que tiende a cero asintóticamente si \(n\lambda^5 \to 0\), por lo que podemos expresar el siguiente nuevo enunciado:
Proposición. Bajo las mismas condiciones del enunciado anterior, con \(f \in C^2[0, 1]\) y \(n\lambda^5 \to 0\), entonces:
\[ \frac{f_\lambda(x) - f(x)}{\sqrt{\text{var}[f_\lambda(x)]}} \xrightarrow{d} N[0, 1] \]
Llegamos entonces a un primer intervalo de confianza para \(f(x)\):
\[ f_\lambda(x) \pm Z_{\alpha/2}\hat{\sigma}\sqrt{\sum_{i=1}^n w_{in}^2(x, \lambda)} \quad (3.26) \]
Para el caso del estimador de Gasser-Müller tenemos que:
\[ w_{in}(x, \lambda) = \lambda^{-1} \int_{s_{i-1}}^{s_i} K(\lambda^{-1}(x-s))ds, \quad i = 1, \ldots, n \]
Y para el caso del estimador de Nadaraya-Watson,
\[ w_{in}(x, \lambda) = \frac{K(\lambda^{-1}(x - x_i))}{\sum_{j=1}^n K(\lambda^{-1}(x - x_j))}, \quad i = 1, \ldots, n \]
En el intervalo (3.26), \(Z_{\alpha/2}\) es el percentil \(100(1-\alpha)\) de la normal estándar y \(\hat{\sigma}\) es la raíz cuadrada de una estimación \(\hat{\sigma}^2\) de \(\sigma^2\) basada en alguno de los estimadores de la varianza en el modelo (2.1) presentados en la sección 2.2.1.
La principal limitación del intervalo (3.26) es que se construye a partir de un estimador \(f_\lambda\) de \(f\) que se sabe que no es insesgado. Vale recordar que hemos elegido \(\lambda\) de tal manera que se minimice el riesgo, que es una combinación del sesgo y de la varianza. Una opción sería intentar construir un intervalo de confianza que considere el sesgo.
Eubank (1999) propone una solución que implica estimar \(f\) con una parámetro de suavización, digamos \(\lambda_1\), estimar luego \(f''\) usando borde en los extremos. El intervalo tomaría la siguiente forma:
\[ f_{\lambda_1}(x) - \frac{\lambda^2}{2} M_2 \hat{f}_{\lambda_2}''(x) \pm Z_{\alpha/2}\hat{\sigma}\sqrt{\sum_{i=1}^n w_{in}^2(x, \lambda)} \quad (3.27) \]
En esta propuesta se debe además adecuar el cálculo de los pesos \(w_{in}(x, \lambda)\) para considerar la forma estimación de \(f\) y tiene como dificultades la elección de más de un parámetro de suavización y el uso de funciones kernel particulares en los extremos. Por otra parte, el intervalo puede verse como si se centrara en una corrección del estimador para lograr que sea insesgado.
Bowman & Azzalini (1997) prefiere utilizar “bandas de variabilidad” que se construyen usando intervalos de confianza calculados para valores de \(x\), alrededor de \(E[f_\lambda(x)]\) usando un estimador de \(\sigma^2\).
La figura 3.7 muestra los intervalos de confianza construidos usando un estimador kernel y las bandas de confianza de Bowman-Azzalini. Ambos ajustes tienen aproximadamente el mismo número de grados de libertad equivalente
# Datos: hora y CO
hora <- 0:23
CO <- c(2.24, 1.53, 1.80, 0.87, 1.21, 1.02, 1.39, 1.90, 2.49, 3.33, 2.82, 3.34,
3.22, 2.94, 2.13, 2.88, 1.95, 1.51, 1.54, 1.60, 1.22, 1.37, 1.96, 1.66)
# Cargar la librería sm (si no está instalada, instálala con install.packages("sm"))
library(sm)
library(ggplot2)
# Establecer el marco para dos gráficos lado a lado
par(mfrow = c(1, 2))
# Estimación kernel con intervalos de confianza al 95%
#sm.regression(x = hora, y = CO, h = 3, model = "none",
# display = "lines", ci = TRUE,
# lty = c(1, 2, 2), col = c("black", "red", "red"),
# xlab = "Hora", ylab = "CO",
# main = "Intervalo de Confianza (95%)")
# Bandas de variabilidad al 95%
#sm.regression(x = hora, y = CO, h = 3, model = "none",
# display = "lines", ci = TRUE, ci.type = "variability",
# lty = c(1, 3, 3), col = c("black", "blue", "blue"),
# xlab = "Hora", ylab = "CO",
# main = "Bandas de Variabilidad (95%)")
sm.regression(x = hora, y = CO, h = 3, display = "se", col = c("black", "red"),
xlab = "Hora", ylab = "CO")
datos <- data.frame(hora = hora, CO = CO)
ggplot(datos, aes(hora, CO)) +
geom_point() +
geom_smooth(method = "loess", se = TRUE, level = 0.95) +
labs(title = "Regresión no paramétrica con bandas de confianza")
# Datos: hora y CO
hora <- 0:23
CO <- c(2.24, 1.53, 1.80, 0.87, 1.21, 1.02, 1.39, 1.90, 2.49, 3.33, 2.82, 3.34,
3.22, 2.94, 2.13, 2.88, 1.95, 1.51, 1.54, 1.60, 1.22, 1.37, 1.96, 1.66)
# Cargar la librería sm (si no está instalada, instálala con install.packages("sm"))
library(sm)
# Establecer el marco para dos gráficos lado a lado
par(mfrow = c(1, 2))
sm.regression(x = hora, y = CO, h = 3, display = "se", col = c("black", "red"),
xlab = "Hora", ylab = "CO")
El término spline proviene del nombre de una herramienta utilizada por ingenieros y arquitectos navales para trazar curvas suaves sobre planos de navios, bajo este concepto de crear curvas suaves que pasan a través de un conjunto de puntos algunos matemáticos comenzaron a trabajar en ello, para poder representar curvas de manera más flexible y precisa que con las técnicas anteriores, como las interpolaciones lineales..
Interpolación y suavización spline
La interpolación es una técnica matemática que nos permite estimar valores desconocidos en función de datos conocidos. Básicamente, consiste en construir una función o curva que pase a través de puntos específicos para luego calcular valores intermedios con base en esa función.
Existe un diferente metodos de interpolación, pero para nuestro estudio hablaremos de la interpolación polinomica, en donde, si se dispone de un conjunto \(A\) de \(k\) puntos en \(\mathbb{R}^2\) entonces se podrá ajustar a los puntos exactamente un polinomio de grado \(k-1\). Es decir, con dos puntos se puede ajustar una recta y con cuatro puntos una cúbica. En esencia si se cuenta con \(K\) puntos entonces el polonomio que se quiere interpolar tendra la forma
\[\pi_{k-1}=a_0+a_1x+a_3x^2+\dots+a_kx^{k-1}.\] Para determinar completamente el polinomio \(\pi_{k-1}\) basta con encontrar los coeficientes \(a_0,a_1,\dots,a_k\) para lo cual se utilizan las coordenadas \((x_i,y_i),i=1,2,\dots,k\) de los k puntos de A. Para esto se requieren \(k\) ecuaciónes, una para cada coeficiente, las que se construyen usando los \(k\) puntos de \(A\).
Por ejemplo, supongamos que tenemos
\[A=\{(1,2);(2,5);(3,4);(4,7)\}\] Y sea \(\pi_3=f(x)=a_0+a_1x+a_2x^2\) el polinomio cuya gráfica pasa por estos cuatro puntos. Para identificar el polinomio \(\pi_3\) usando interpolación, basta con encontrar los valores de las constantes \(a_0,a_1,a_2\).
Al resolver obtenemos que los valores son: \((-13,23.66,-10)\) y el gráfico del polinomio aparece a continuación.
#> (Intercept) poly(x, degree = 3, raw = TRUE)1
#> -13.00000 23.66667
#> poly(x, degree = 3, raw = TRUE)2
#> -10.00000# Secuencia de puntos para graficar el polinomio
x_vals <- seq(1, 4, by = 0.01)
y_vals <- predict(modelo, newdata = data.frame(x = x_vals))
# Graficar puntos originales y curva interpolada
plot(x, y, col = "red", pch = 19, main = "Interpolación Polinómica de Grado 3", xlab = "x", ylab = "y")
lines(x_vals, y_vals, col = "blue", lwd = 2)
Una propiedad llamativa de un polinomio cúbico es que tiene un único punto de inflexión. Si se desea construir una curva que tiene más puntos de inflexión, una solución posible es resolver un sistema más grande de ecuaciones lineales. Otra solución posible en este caso sería construir la interpolación uniendo varios polinomios cúbicos.
Para lograrlo se escogen varios puntos llamados comúnmente . Sean \((p_1, p_2, \ldots, p_n)\) los puntos de control y \((f_1, f_2, \ldots, f_{n-1})\) las cúbicas que se usarán para interpolar \(f\).
El propósito es ajustar una cúbica entre cada par de puntos de control \((p_i, p_{i+1})\). Si se denota $ f_i(x) = a_{i0} + a_{i1}x + a_{i2}x^2 + a_{i3}x^3 $ a la cúbica que une estos dos puntos de control.
Para construir estas ecuaciones es necesario imponer algunas condiciones mínimas a las cúbicas \(f_i\) que se unirán entre sí para construir el polinomio \(f\).
Estimación
Un problema igualmente interesante desde el punto de vista numérico pero más interesante desde el punto de vista estadístico es aquel en el cual se sabe que los datos disponibles son mediciones sujetas a error y que aunque se cree que la relación entre \(X\) y \(Y\) es funcional, no necesariamente los puntos observados están sobre la curva que relaciona las dos variables, como se ve en la siguiente grafica
Entonces, el propósito de la estimación es proponer una función \(\hat{f}\) que permita aproximar \(f\).
En realidad en ambos casos, en la interpolación y en la estimación, se obtienen curvas aproximadas. La diferencia está en que, por una parte, no es posible utilizar interpolación con los puntos y que a través de los métodos de estimación es posible asociar a las curvas construidas una medida probabilística de qué tan buena será la aproximación alcanzada.
Una posible solución estadística es la regresión polinómica de la sección @ref(subsec-polinomica), en la que se asume que la función subyacente \(f\) es una función continua que se aproxima mediante un polinomio de Taylor y que las \(n\) observaciones disponibles provienen de tal función \(f\), de tal manera que los valores \(y_i\) se separan de los de \(f(x_i)\) por una cantidad aleatoria \(\varepsilon_i\).
Estimación spline
Una función spline es una función suave a tramos que se utiliza para aproximar a a un conjunto de datos atraves de splines, estos unen varios polinomios de bajo grado en distintos intervalos, asegurando que la transición entre ellos sea suave.
Supongamos que tenemos un conjunto de números reales \(x_1, \ldots, x_n\) en un intervalo \([a, b]\), tales que \[ a < x_1 < x_2 < \ldots < x_n < b. \]
Una función \(s\) definida en \([a, b]\) es un spline cúbico si cumple las siguientes dos condiciones:
A los puntos \(x_i\) los llamaremos nodos.
Definición Un spline cúbico en \([a, b]\) se llama un spline cúbico natural si se cumple que las dos primeras derivadas de \(s\) son iguales a cero en los puntos extremos \(a\) y \(b\). Estas condiciones se conocen como condiciones de acotamiento natural, e implican que \(s\) es lineal en los dos intervalos extremos \((a, x_1)\) y \((x_n, b)\).
En la estimación usando series de Fourier acordamos que la función de regresión está en \(L_2[0, 1]\). Añadiremos ahora la condición de que las dos primeras derivadas de \(f\), que denotaremos \(f'\) y \(f''\), también pertenezcan a \(L_2[0, 1]\). A este espacio de funciones lo llamaremos el espacio de Sobolev de orden 2 y lo denotaremos por \[ W_2^2[0, 1]. \]
Una medida natural de suavidad asociada con una función \(f \in W_2^2[0, 1]\) es \[ \int_0^1 \left( f''(x) \right)^2 \, dx, \] mientras que una medida de bondad de ajuste de los datos al modelo es la suma de cuadrados del error: \[ \frac{1}{n} \sum_{i=1}^{n} \left( y_i - f(x_i) \right)^2. \] Esto implica que una medida de la calidad de un estimador de \(f\) podría basarse en la siguiente suma convexa: \[ (1 - q)\frac{1}{n} \sum_{i=1}^{n} \left( y_i - f(x_i) \right)^2 + q \int_0^1 \left( f''(x) \right)^2 dx, \] con \(0 < q < 1\).
Si hacemos \(\lambda = \dfrac{q}{1 - q}\), la elección del estimador de \(f\) es equivalente a elegir \(f_\lambda\) que minimice la siguiente expresión:
\[\begin{equation} \frac{1}{n} \sum_{i=1}^{n} \left( y_i - f(x_i) \right)^2 + \lambda \int_0^1 \left( f''(x) \right)^2 dx, (\#eq:spline) \end{equation}\] sobre todas las funciones \(f \in W_2^2[0, 1]\). A este estimador \(f_\lambda\) lo llamaremos un estimador spline de \(f\).
De la expresión @ref(eq:spline) se sigue que:
Eubank en su libro “Nonparametric Regression and Spline Smoothing” demostro que la solución a este problema de optimización es única, y corresponde al estimador: \[ f_\lambda = \sum_{i=1}^{n} \beta_{\lambda i} f_i, \] donde los \(\beta_{\lambda i}\) son coeficientes reales, y \(\{f_i\}\) es un conjunto de funciones base..
Donde \(\boldsymbol{\beta}_\lambda = (\beta_{\lambda 1}, \beta_{\lambda 2}, \ldots, \beta_{\lambda n})^\top\) es la única solución con respecto a \[ \mathbf{c} = (c_1, c_2, \ldots, c_n)^\top \] del sistema de ecuaciones: \[\begin{equation} (X^\top X + n\lambda \, \Omega) \, \mathbf{c} = X^\top \mathbf{y} \end{equation}\] donde: \[ X = \{f_j(x_i)\}_{i,j=1,2,\ldots,n}, \quad \mathbf{y} = (y_1, y_2, \ldots, y_n)^\top, \] y \[ \Omega = \left[ \int_0^1 f_i''(x) f_j''(x) \, dx \right]_{i,j=1,2,\ldots,n}. \]
Las funciones \(\{f_j\}_{j=1}^{n}\) forman una base del conjunto de splines cúbicos naturales.
Si usamos la base de splines cúbicos naturales, entonces se tendría que el vector de valores estimados es \[ \mathbf{f}_\lambda = \begin{pmatrix} f_\lambda(x_1), f_\lambda(x_2), \ldots, f_\lambda(x_n) \end{pmatrix}^T = S_\lambda \mathbf{y}, \] donde \[\begin{equation} S_\lambda = X (X^T X + n\lambda \Omega)^{-1} X^T. (\#eq:spline2) \end{equation}\]
La elección del parámetro de suavización \(\lambda\) se hace usualmente con el estimador de validación cruzada generalizada GCV, usando la matriz \(S_\lambda\).
En R existe la función smooth.spline del paquete stats que selecciona \(\lambda\) y realiza la estimación.
#datos
hora_reescalada <- (0:23)/23
co <- c(2.24, 1.53, 1.80, 0.87, 1.21, 1.02, 1.39, 1.90, 2.49, 3.33,
2.82, 3.34, 3.22, 2.94, 2.13, 2.88, 1.95, 1.51, 1.54, 1.60,
1.22, 1.37, 1.96, 1.66)
# Ajuste sin spar
spline <- smooth.spline(hora_reescalada, co)
print(spline)#> Call:
#> smooth.spline(x = hora_reescalada, y = co)
#>
#> Smoothing Parameter spar= 0.5058871 lambda= 0.0002605305 (14 iterations)
#> Equivalent Degrees of Freedom (Df): 7.306994
#> Penalized Criterion (RSS): 1.577585
#> GCV: 0.1358735
Aunque tambien la función te permite modificar el parámetro spar en la función smooth.spline() de R controla el grado de suavizado del ajuste del spline cúbico suavizado, este valor esta entre 0 y 1, donde 1 es la maxima suavización.
#Datos
hora_reescalada <- (0:23)/23
co <- c(2.24, 1.53, 1.80, 0.87, 1.21, 1.02, 1.39, 1.90, 2.49, 3.33,
2.82, 3.34, 3.22, 2.94, 2.13, 2.88, 1.95, 1.51, 1.54, 1.60,
1.22, 1.37, 1.96, 1.66)
# Ajuste sin spar
spline_1 <- smooth.spline(hora_reescalada, co,spar = 0.95)
spline_2 <- smooth.spline(hora_reescalada, co,spar = 0.4)
# Gráfico
plot(hora_reescalada, co, col="blue", pch=19)
lines(predict(spline_1, seq(0,1,0.01)), col="red", lwd=2)
lines(predict(spline_2, seq(0,1,0.01)), col="purple", lwd=2)
legend("topright", legend = c("spar = 0.95", "spar = 0.4"),
col = c("red", "purple"), lwd = 1, bty = "n", cex = 1)
Estimación spline por minimos cuadrados
Las propiedades asintóticas de los estimadores spline que hemos discutido hasta aquí se derivan del hecho de que estos estimadores pueden representarse como estimadores kernel, de tal manera que sus propiedades asintóticas son similares.
Existe, sin embargo, otro posible estimador spline que no puede representarse como un estimador kernel. Esta nueva forma de estimación se apoya en el modelo de regresión polinómica discutido. En ese caso, representamos el modelo @ref(eq:Reglin) de la siguiente manera:
\[\begin{equation} f(x_i) = \sum_{j=0}^{\lambda} \theta_j f_j(x_i) + r_\lambda(x_i), \quad i = 1, \ldots, n (\#eq:splinec) \end{equation}\]
donde \(\lambda\) cumple el papel de parámetro de suavización, y consiste en el número de términos que se conservan en la sumatoria.
Si modificamos la notación y llamamos \(m\) a lo que en la expresión @ref(eq:splinec) hemos denotado como \(\lambda\), una de las representaciones posibles para \(r_m(x)\) sería:
\[\begin{equation} r_m(x) = \frac{1}{m - 1} \int_0^1 f^{(m)}(\xi) (x - \xi)^{m-1}_+ \, d\xi \end{equation}\]
Si \(r(x_1), \ldots, r(x_n)\) son muy pequeños, entonces el modelo de regresión polinómica sería adecuado para un conjunto particular de datos que tenga estas características. Pero si este no es el caso, una posible generalización es la , que se apoya en la siguiente aproximación a la integral:
\[\begin{equation} r_m(x) \approx \sum_{j=1}^{k} \delta_j (x - \xi_j)^{m-1}_+ \end{equation}\]
para un conjunto de constantes \(\delta_1, \ldots, \delta_k\) y un conjunto de puntos \(0 < \xi_1 < \ldots < \xi_k < 1\), de tal manera que una nueva aproximación a la función de regresión lucirá ahora así:
\[\begin{equation} f(x) = \sum_{j=1}^{\lambda} \theta_j x^{j-1} + \sum_{j=1}^{k} \delta_j (x - \xi_j)^{m-1}_+ \end{equation}\]
Sea \(\{\xi_1, \ldots, \xi_k\}\) un conjunto dado de puntos (nodos). Una posible solución para estimar \(f\) sería estimar por mínimos cuadrados los \(m\) coeficientes \(\{\theta_j\}_{j=1}^{m}\) y los \(k\) coeficientes \(\{\delta_j\}_{j=1}^{k}\).
Sea entonces el vector de parámetros:
\[ \boldsymbol{\beta} = (\theta_1, \ldots, \theta_m, \delta_1, \ldots, \delta_k)^\top \]
El estimador spline por mínimos cuadrados de \(f\) será:
\[\begin{equation} f_\lambda(x) = \sum_{j=1}^{m+k} \beta_{\lambda j} f_j(x) \end{equation}\]
Un problema aún en curso es la elección de \(m\) y de \(k\). A los \(k\) valores de \(\xi\) se les llama nodos, y la solución más sencilla es hacer uso de inspección visual de los datos, cuando sea posible.
Si \(f\) cambia muy rápidamente, entonces deberemos usar más nodos. La elección de \(m\) depende de si se usará un spline lineal, cuadrático, cúbico, etc., y los \(k\) nodos se ubicarán consecuentemente.
En R existe la funcíon bs para ajustar splines con mínimos cuadrados ordinarios, la función genera una base de B-splines (Basis Splines) para ajustar modelos de splines.
library(splines)
# Datos
x <- (0:23)/23
y <- c(2.24, 1.53, 1.80, 0.87, 1.21, 1.02, 1.39, 1.90,
2.49, 3.33, 2.82, 3.34, 3.22, 2.94, 2.13, 2.88,
1.95, 1.51, 1.54, 1.60, 1.22, 1.37, 1.96, 1.66)
# Definir los nodos
knots <- quantile(x, probs = seq(0.1, 0.9, length.out = 5))
# Crear matriz B-splines cúbicos
X <- bs(x, knots = knots, degree = 3, intercept = TRUE)
# Ajustar modelo por mínimos cuadrados
fit <- lm(y ~ X - 1) # EL intercepto ya esta incluido en bs()
# Crear puntos para la curva ajustada
x_new <- seq(0, 1, length.out = 100)
X_new <- predict(bs(x, knots = knots, degree = 3, intercept = TRUE), x_new)
y_pred <- X_new %*% coef(fit)
# Graficar resultados
plot(x, y, pch = 19, col = "blue", main = "B-splines")
lines(x_new, y_pred, col = "red", lwd = 2)
points(knots, predict(fit, newdata = list(X = predict(bs(x, knots = knots, degree = 3, intercept = TRUE), knots))),
pch = "x", col = "green", cex = 1.5)
legend("topright", legend = c("Ajuste", "Nodos"),
col = c("red", "green"), pch = c(NA, "x"), lwd = c(2, NA))
Los Modelos Aditivos Generalizados (GAMs) extienden los modelos lineales permitiendo relaciones no lineales entre predictores y la variable respuesta, manteniendo una estructura aditiva.
Para una variable respuesta \(Y\) con distribución en la familia exponencial y predictores \(X_1,\ldots,X_p\), el GAM se define como:
\[ g\left(\mathbb{E}[Y|X]\right) = \beta_0 + \sum_{j=1}^p f_j(X_j) \]
donde:
\[ \begin{aligned} \text{Función de enlace:} \quad & \eta = g(\mu) \\ \text{Función de varianza:} \quad & \text{Var}(Y) = V(\mu)\phi \\ \text{Penalización:} \quad & \int \left[f_j''(x)\right]^2 dx \end{aligned} \]
Se maximiza la log-verosimilitud penalizada, para mas detalles Clemente (2022)
\[ \ell_p(\boldsymbol{\theta}) = \ell(\boldsymbol{\theta}) - \frac{1}{2}\sum_{j=1}^p \lambda_j \int \left[f_j''(x)\right]^2 dx \]
donde \(\boldsymbol{\theta}\) incluye \(\beta_0\) y los parámetros de las funciones suavizadas \(f_j\).
Familias Comunes
| Familia | Función de enlace | Aplicación |
|---|---|---|
| Gaussiana | Identidad | Datos continuos |
| Gamma | Log | Datos positivos |
| Poisson | Log | Conteos |
| Binomial | Logit | Proporciones |
A modo de compración:
| Característica | GLM | GAM |
|---|---|---|
| Forma funcional | Lineal (\(\beta_j X_j\)) | No lineal (\(f_j(X_j)\)) |
| Término error | \(\varepsilon\) explícito | Implícito en distribución |
Para los datos en la sección @(Modelos paramétricos)
library(ggplot2)
library(mgcv)
set.seed(123)
x <- seq(1, 10, length.out = 100)
y <- exp(0.3 * x + rnorm(100, sd = 0.2))
datos <- data.frame(x = x, y = y)
modelo1 <- gam(y ~ s(x), data = datos)
modelo2 <- gam(y ~ s(x, k = 15), family = Gamma(link = "log"), method = "REML", data = datos)
AIC(modelo1, modelo2)#> df AIC
#> modelo1 6.022066 404.9719
#> modelo2 3.003534 282.5270summary(modelo1)#>
#> Family: gaussian
#> Link function: identity
#>
#> Formula:
#> y ~ s(x)
#>
#> Parametric coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 7.3333 0.1771 41.41 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Approximate significance of smooth terms:
#> edf Ref.df F p-value
#> s(x) 4.022 4.965 208.1 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> R-sq.(adj) = 0.913 Deviance explained = 91.6%
#> GCV = 3.3017 Scale est. = 3.1359 n = 100summary(modelo2)#>
#> Family: Gamma
#> Link function: log
#>
#> Formula:
#> y ~ s(x, k = 15)
#>
#> Parametric coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 1.68453 0.01851 91.01 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Approximate significance of smooth terms:
#> edf Ref.df F p-value
#> s(x) 1.002 1.004 1868 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> R-sq.(adj) = 0.912 Deviance explained = 94.9%
#> -REML = 144.44 Scale est. = 0.034258 n = 100
#>
#> Method: GCV Optimizer: magic
#> Smoothing parameter selection converged after 6 iterations.
#> The RMS GCV score gradient at convergence was 3.073791e-07 .
#> The Hessian was positive definite.
#> Model rank = 10 / 10
#>
#> Basis dimension (k) checking results. Low p-value (k-index<1) may
#> indicate that k is too low, especially if edf is close to k'.
#>
#> k' edf k-index p-value
#> s(x) 9.00 4.02 1.03 0.68
#>
#> Method: REML Optimizer: outer newton
#> full convergence after 9 iterations.
#> Gradient range [-7.231766e-05,0.0007281679]
#> (score 144.4377 & scale 0.03425775).
#> Hessian positive definite, eigenvalue range [7.344579e-05,49.55541].
#> Model rank = 15 / 15
#>
#> Basis dimension (k) checking results. Low p-value (k-index<1) may
#> indicate that k is too low, especially if edf is close to k'.
#>
#> k' edf k-index p-value
#> s(x) 14 1 1.03 0.59
Definición: EL espacio \(L_2[0,1]\) se define como el espacio funcional de todas las funciones cuadrado integrables en [0, 1].
Sean \(f_1\) y \(f_2\) dos funciones en \(L_2[0, 1]\). Entonces, la norma de \(f_i\) se define como:
\[ \|f_i\| = \left( \int_0^1 f_i^2(x) \, dx \right)^{1/2}, \quad i = 1, 2 \] mientras que el producto interno de \(f_1\) y \(f_2\) es:
\[ <f_1, f_2 > = \int_0^1 f_1(x) f_2(x) \, dx \]
En un espacio de dimensión finita, por ejemplo en \(\mathbb{R}^p\), siempre es posible representar cualquier elemento del espacio usando una combinación lineal de elementos en una base del espacio.
Pero el espacio \(L_2[0, 1]\) es un espacio de dimensión infinita, por lo que las siguientes definiciones se requieren para poder concluir que, en efecto, es posible representar cualquier elemento de \(L_2[0, 1]\) a partir de un conjunto de elementos que formen una base de \(L_2[0, 1]\).
Dos funciones \(f_1, f_2 \in L_2[0, 1]\) se dicen ortogonales si \(\langle f_1, f_2 \rangle = 0\).
La ortogonalidad de \(f_1\) y \(f_2\) se denota \(f_1 \perp f_2\).
Una sucesión de funciones \(\{f_j\}_{j=1}^{\infty} \subset L_2[0, 1]\) se dice ortonormal
si las \(f_j\) son ortogonales por parejas y \(\|f_j\| = 1\) para todo \(j\).
Una sucesión de funciones \(\{f_j\}_{j=1}^{\infty} \subset L_2[0, 1]\) se dice que es una
sucesión ortonormal completa (CONS, por sus siglas en inglés) si
\(f \perp f_j\) para todo \(j\) implica que \(f = 0\).
Las definiciones de ortogonalidad y ortonormalidad en \(L_2[0, 1]\) son similares a las mismas definiciones en \(\mathbb{R}^p\), mientras que la definición de una CONS implica que la única función ortogonal a todas las funciones \(f_j\) es la función cero.
Coeficientes de Fourier
La siguinete proposición nos permite concluir que, cualquier CONS provee una base de \(L_2[0, 1]\). Por lo que una función \(f\) de \(L_2[0, 1]\) puede representarse como una combinación lineal de una colección de funciones \(\{f_j\}_{j=1}^\infty\).
Proposicion Sea \(\{f_j\}_{j=1}^{\infty}\) una CONS de \(L_2[0,1]\) y sea \(f\) una función de \(L_2[0,1]\). Si definimos \[\begin{equation} f_\lambda(x)=\sum_{i=1}^nK(x, x_i; \lambda)y_i (\#eq:CofFour) \end{equation}\] entonces \(\sum_{j=1}^{\lambda} \beta_j f_j\) es la mejor aproximación a \(f\) en el sentido de que \[ \left\| f - \sum_{j=1}^{\lambda} \beta_j f_j \right\| \leq \left\| f - \sum_{j=1}^{\lambda} b_j f_j \right\| \] para todo \(b = (b_1, b_2, \dots, b_\lambda)^T \in \mathbb{R}^\lambda\). Más aún, cuando \(\lambda \to \infty\), \[ \left\| f - \sum_{j=1}^{\lambda} \beta_j f_j \right\| \to 0. \]
Llamaremos coeficientes generalizados de Fourier a los coeficientes \(\beta_j\) de la expresión @ref(eq:CofFour).
Finalmente, se sabe que los coeficientes generalizados de Fourier satisfacen la relación de Parseval. \[\sum_{j=1}^\infty \beta_j^2=\|f\|.\]
Definicón Dada una secuencia de datos \(\{y_i\}_{i=1}^n\) una diferencia de orden \(r\) con pesos \(\{d_k\}_{k=-m_1}^{m_2}\), se calcula como: \[ \Delta_l=\sum_{k = -m_1}^{m_2} d_k y_{k + l}. \] donde los \(d_k\) satisfacen: \(\sum_{k=-m_1}^{m_2} d_k = 0\) (elimina tenedencias de bajo grado) y \(\sum_{k=-m_1}^{m_2} d_k^2 = 1\) (mantiene la escala de la varianza).
Dado un conjunto de nodos: \[ a = x_0 < x_1 < x_2 < \cdots < x_n = b, \] un spline cúbico \(S(x)\) está definido como: \[ S(x) = \begin{cases} P_1(x), & x \in [x_0, x_1] \\ P_2(x), & x \in [x_1, x_2] \\ \vdots \\ P_n(x), & x \in [x_{n-1}, x_n] \end{cases} \] donde cada \(P_i(x)\) es un polinomio cúbico de la forma: \[ P_i(x) = a_i + b_i(x - x_{i-1}) + c_i(x - x_{i-1})^2 + d_i(x - x_{i-1})^3 \]
Equipo 4:
Oswaldo Bueno
Alejandra Borrego
Edruan Ramírez
Equipo 5:
Juan Pablo Rodriguez
Luz María Salazar
Héctor Soto