El término ‘spline’ se refiere a la herramienta de un artesano, una tira fina y flexible de madera o metal, utilizada para trazar curvas suaves. Se aplicaban varios pesos en varias posiciones para que la tira se doblara según su número y posición. Ésta se obligaba a pasar por un conjunto de puntos fijos: clavijas de metal, las costillas de un barco, etc. En una superficie plana solían ser pesas con un gancho incorporado y, por tanto, fáciles de manipular. La forma del material doblado adoptaría naturalmente la forma de una curva spline. Del mismo modo, las splines se utilizan en estadística para reproducir matemáticamente las formas flexibles. Los nudos se colocan en varios lugares dentro de la gama de datos, para identificar los puntos en los que las piezas funcionales adyacentes se unen entre sí. En lugar de estrías de metal o madera, se eligen piezas funcionales suaves (normalmente polinomios de bajo orden) para ajustar los datos entre dos nudos consecutivos. El tipo de polinomio y el número y la colocación de los nudos es lo que define entonces el tipo de spline.
Ejemplo motivador
Con la introducción de los modelos aditivos generalizados (GAM) en 1986, el uso del modelado spline se ha convertido en una herramienta establecida en el análisis de regresión estadística. Para ilustrar esto, consideremos los datos de un conjunto de 892 mujeres menores de 50 años recogidos en tres pueblos de África Occidental (datos disponibles en el archivo adicional 1: Apéndice). Queremos explorar la relación entre la edad (en años) y una medida cruda de la grasa corporal, que es el grosor del pliegue cutáneo del tríceps. La figura 1 muestra la relación entre la edad y el grosor del pliegue cutáneo del tríceps medido en escala logarítmica. La figura 1 muestra la relación entre la edad y el grosor del pliegue cutáneo del tríceps medido en escala logarítmica. La línea discontinua representa un ajuste lineal simple, la línea sólida un ajuste utilizando polinomios flexibles de tercer grado
Un modelo de regresión simple de la forma yi=β0+β1xi+ε,i=1,…,n, difícilmente daría una aproximación al patrón observado, ya que es obvio que la relación no es lineal. El modelo puede ampliarse para dar cabida a los efectos no lineales utilizando algunos polinomios. Entonces, los efectos no lineales podrían ser modelados por un polinomio de grado 3 dado por:
donde u es una función de x llamada función base, definida aquí por:
El modelo de regresión descrito en la Ec. 1 sigue siendo un modelo lineal, a pesar de que proporciona una función no lineal de la variable predictora. El modelo sigue siendo lineal en los coeficientes y puede ajustarse utilizando métodos de mínimos cuadrados ordinarios. La base puede crearse en R utilizando la función poly(x,3) con las entradas x (referida a la variable), y p (referida al grado del polinomio). Esto conduce a un modelo simple y suave univariante de la forma: yi=f(xi)+ε donde f() es alguna función/transformación del predictor. Este modelo puede ajustarse fácilmente en R utilizando: lm(y ∼poly(x,3)). A pesar de la simplicidad, la regresión polinómica tiene varios inconvenientes, el más importante es la no localización. Esto significa que la función ajustada en un valor dado x0 depende de valores de datos alejados de ese punto. Es fácil ver esto en acción ajustando un polinomio a un conjunto de datos y moviendo uno de los puntos de datos cerca del borde derecho hacia arriba o hacia abajo. Como resultado, la función ajustada normalmente cambiará lejos de esa coordenada x.
Considere, en lugar de ajustar un polinomio global, particionar el rango de x en intervalos más pequeños, utilizando un número y posición arbitrarios de puntos, τ, también llamados los nodos. Se puede ajustar un modelo continuo simple a trozos definiendo las funciones: f1(x)=1,f2(x)=x,f3(x)=(x-τ1)+,f4(x)=(x-τ2)+,…, siendo «+» una función definida como: $$
El conjunto de estas funciones conducen a una función compuesta f(x).
Definición de splines
La spline metálica del dibujante puede asumir formas arbitrarias, por ejemplo, la sección transversal de un ala de avión o la espiral de una bomba centrífuga. Para las aplicaciones estadísticas asumiremos curvas de la forma f(X), es decir, un único valor y para cada x. El predictor x puede ser una única variable o múltiples variables. Nuestra discusión se centrará casi por completo en una función univariante con \(X\ en \mathbb {R}\). Definamos un conjunto de nudos τ1<…<τK en el rango de X. Un spline f(X) será una función suave, que satisface ciertas propiedades de diferenciabilidad que se mencionan más adelante, tal que f(X) es un polinomio de grado d. Los splines de madera o metal tienen derivadas continuas de todos los órdenes ya que son un objeto físico. Esto no es cierto para los splines estadísticos. Más bien imponemos un criterio de suavidad según el cual todas las derivadas de orden inferior a d son continuas. Un spline físico es lineal más allá del último nudo y podemos imponer una restricción adicional de que las derivadas de orden 2 o superior sean nulas en los nudos más a la izquierda y a la derecha; los splines con esta restricción adicional se conocen como splines «restringidos» o «naturales». Para obtener curvas más flexibles se puede aumentar el número de nudos o el grado del polinomio. Sin embargo, hay una compensación; aumentar el número de nudos puede sobreajustar los datos y aumentar la varianza, mientras que disminuir el número de nudos puede dar lugar a una función rígida y restrictiva que tiene más sesgo.
Representación mediante funciones base
Supongamos que la función desconocida f está representada por una función spline con secuencia de nudos fija y grado fijo d. Dado que estas últimas funciones forman un espacio vectorial V, es posible escribir f como
donde las Bk son un conjunto de funciones base que definen V y βk son los coeficientes del spline asociados. Con k nudos hay k+1 polinomios de grado d junto con d∗k restricciones, lo que lleva a (d+1)(k+1)-d∗k=d+k+1 parámetros libres; para un spline natural hay k parámetros libres. Dado que βB=(βA)(A-1B)=γB∗ para cualquier matriz no singular A hay un número infinito de conjuntos de bases posibles para el ajuste del spline.
La representación en (2) tiene la ventaja de que la estimación de f se reduce a la estimación de los coeficientes βk. Más concretamente, la expresión en (2) es lineal en el vector de coeficientes β=(β1,…,βK+d+1). Por lo tanto, la estimación de f puede verse como un problema de optimización que es lineal en las variables transformadas B1(X),…,BK+d+1(X), lo que permite el uso de técnicas de estimación bien establecidas para el uso de splines en una amplia gama de modelos de regresión multivariable (generalizada). Es importante destacar que el modelado spline reduce la estimación de las funciones f() a la estimación de un pequeño conjunto de coeficientes de valor real.
Como han señalado varios autores (por ejemplo, la gran flexibilidad del modelado spline tiene el precio de una serie de parámetros de ajuste. Dos de ellos, la elección de las funciones base B y el grado d de los polinomios subyacentes resultan tener poco impacto. Los polinomios cúbicos (d=3) son la norma habitual, ya que dan lugar a curvas que parecen perfectamente suaves para el ojo humano. Si las derivadas de las curvas ajustadas son de interés, un orden más alto es a veces apropiado, pero en general los ajustes para d>3 son efectivamente indistinguibles. Los ajustes con d=1 o d=2 tienen propiedades estadísticas casi idénticas, pero parecerán más irregulares. La elección entre dos conjuntos de bases B y B∗ no cambiará, por definición, las predicciones de un ajuste, por lo que se reduce a cuestiones de conveniencia.
Las dos opciones clave están en el número y el espaciado de los nudos y el uso (o no) de una función de penalización, por ejemplo, la segunda derivada integrada del spline. Cuando no hay penalización, la creación de las variables transformadas puede hacerse por separado y las nuevas variables se incluyen simplemente en un ajuste de modelo estándar; no es necesario modificar el procedimiento de regresión subyacente. Este enfoque suele denominarse splines de regresión; la flexibilidad de la función no lineal resultante es totalmente una función del número de nudos. En cambio, la inclusión de una penalización por suavización requiere la modificación de la rutina de ajuste para poder acomodarla. Ésta debe incluirse en cada función de regresión por separado. Los splines de suavizado resultantes tienen varias propiedades deseables, pero la complejidad añadida de la función de suavizado puede ser una razón para que no se utilicen más a menudo en entornos aplicados.
Aunque se ha llevado a cabo una investigación considerable para explorar las propiedades matemáticas de los diversos enfoques de splines (véase , los estadísticos aplicados y los analistas de datos apenas parecen ser conscientes de estos resultados cuando utilizan el modelado de splines en aplicaciones prácticas. De hecho, muchos de los artículos identificados por nuestra búsqueda en la web no contenían ninguna justificación sobre el fundamento de la elección del método spline utilizado.
Bases spline populares
Existen numerosas opciones para la definición de las funciones base Bk, donde las diversas bases spline difieren con respecto a sus propiedades numéricas . En esta sección, vamos a introducir algunas de las bases spline más populares, a saber, la base de la serie de potencia truncada, la base B-spline y la base spline cardinal.
Serie de potencia truncada y splines cúbicos
La base de la serie de potencia truncada está definida por las funciones de base
Una ventaja de las funciones base anteriores es su fácil interpretación: Partiendo de un polinomio «básico» de grado d definido en (primera línea de la ecuación), se añaden sucesivamente desviaciones del polinomio básico a la función spline a la derecha de cada uno de los K nudos (segunda línea). Un spline de base de potencia truncada es d-1 veces diferenciable en los nudos y tiene d+K grados de libertad. Es relativamente fácil para el usuario crear una serie de potencia truncada en R. Dejemos que x represente algunas observaciones en , entonces se puede crear una base de potencia truncada de grado d=3 con 5 nudos igualmente espaciados a lo largo del rango de x utilizando el Código 1 en el Archivo adicional 1: Apéndice (Fig. 2).
Una característica de las series de potencia truncadas es que los soportes de las funciones no son locales, con algunos de los Bk definidos en todo el rango de datos . Esto podría conducir a altas correlaciones entre algunas splines de base, lo que implica inestabilidades numéricas en la estimación de splines. Para la base de serie de potencia truncada, se da un ejemplo en , Capítulo 5.
Los splines cúbicos se crean utilizando un polinomio cúbico en un intervalo entre dos nudos sucesivos. El spline tiene cuatro parámetros en cada una de las regiones K+1 menos tres restricciones para cada nudo, lo que resulta en un K+4 grados de libertad.
Una función spline cúbica, con tres nudos (τ1,τ2,τ3) tendrá 7 grados de libertad. Utilizando la representación dada en la Ecuación 2, la función se puede escribir como:
B-splines
La base B-spline es una base spline comúnmente utilizada que se basa en una parametrización especial de un spline cúbico. La base B-spline , se basa en la secuencia de nudos
donde los conjuntos ξd+2 := τ1,…,ξd+K+1:=τK y ξd+1:=a,ξd+K+2:=b se denominan «nudos interiores» y «nudos límites», respectivamente. La elección de los nudos adicionales ξ1,…,ξd y ξd+K+3,…,ξ2d+K+2 es esencialmente arbitraria. Una estrategia común es establecerlos igual a los nudos límite. Alternativamente, si los nudos interiores y los nudos límite ξd+1<…<ξd+K+2 se eligen para que sean equidistantes, es decir, ξk+1-ξk=δ ∀k∈{d+1,…,d+K+1}, los nudos límite pueden situarse en ξd+1-δ,…,ξd+1-d-δ y ξd+K+2+δ,…,ξd+K+2+d-δ.
Para d>0, las funciones de base B-spline de grado d (denotadas por \(B_{k}^{d}(x)\N) se definen mediante la fórmula recursivaFootnote 1
donde
y \(B_{k}^{0}(x) \equiv 0\) si ξk=ξk+1. Las B-splines tienen la ventaja de que las funciones base tienen soporte local. Más concretamente, son mayores que cero en los intervalos abarcados por d+2 nudos y cero en el resto. Esta propiedad resulta en una alta estabilidad numérica, y también en un algoritmo eficiente para la construcción de las funciones base, ver para más detalles.
Splines naturales cúbicos y cardinales
Un spline polinómico como un cúbico o un B-spline, puede ser errático en los límites de los datos. Para solucionar este problema, los splines naturales son splines cúbicos que tienen la restricción adicional de que son lineales en las colas de los nudos de los límites (-∞,a],, Capítulo 4.
Además de las series de potencia truncadas splines naturales, B-spline y bases de splines cardinales, existen otras bases menos populares. Para una visión general, nos remitimos a los libros de.
Splines penalizados
Los splines presentados hasta ahora suelen denominarse splines de regresión. Además de la elección de la base del spline (B-spline, serie de potencia truncada, etc.), hay que elegir el número de nudos y las posiciones de los mismos. Obviamente, estos parámetros de ajuste pueden tener un impacto importante en la forma estimada de una función spline: Un gran número de nudos implica una gran flexibilidad, pero también puede dar lugar a un ajuste excesivo de los datos en cuestión. Por el contrario, un número reducido de nudos puede dar lugar a una estimación «excesivamente suave» que es propensa a un sesgo de infraajuste (véase ).
Un enfoque popular para facilitar la elección de las posiciones de los nudos en el modelado de splines es el uso de splines penalizados. Dada una muestra i.i.d. de datos (x1,y1),…(xn,yn), un spline penalizado es la solución al problema
donde lβ denota la log-verosimilitud (o, en el caso de la regresión de Cox, la log-verosimilitud parcial) y Jr es una penalización de rugosidad que se hace pequeña si la función spline es «suave». En general, los splines penalizados se basan en la idea de que la función desconocida f se modela mediante un spline con un gran número de nudos, lo que permite un alto grado de flexibilidad. Por otro lado, una estimación aproximada del spline que tenga un valor alto de lβ y se acerque a los valores de los datos da lugar a un valor grande de Jβ. La maximización de esta función implica, por tanto, un compromiso entre la suavidad y el ajuste del modelo que es controlado por el parámetro de ajuste λ≥0.
Un caso especial es el problema de mínimos cuadrados penalizados
en la regresión gaussiana. La penalización \ {J_{beta } \N, \Nint _{a}^{b} \left (\partial ^{2} f / \partial x^{2}\right)^{2} dx\) expresa la «suavidad» de una función spline en términos de la segunda derivada de f. Para un λ dado, se puede demostrar que la solución es un spline cúbico natural con la secuencia de nudos x(1)<…<x(n), es decir, las posiciones de los nudos no tienen que ser elegidas, sino que están dadas «naturalmente» por los valores de datos únicos ordenados de X. En la literatura, este tipo de spline se denomina spline de suavizado. Cabe destacar que se puede demostrar que un spline de suavizado interpola los datos si λ=0, mientras que λ=∞ implica una función lineal. Obsérvese que los splines de suavizado son un caso especial de la clase más general de splines de placa fina , que permiten una extensión del criterio de la Ec. (3) a xi de mayor dimensión (véase , Sección 4.15], y para más detalles).
Una propiedad conveniente de los splines de suavizado es que la penalización Jβ puede escribirse como β⊤Ωβ con una matriz de penalización Ω convenientemente definida. Por lo tanto, la solución a (3) viene dada por la estimación de mínimos cuadrados penalizados
donde B es una matriz de dimensión n×n que contiene las funciones base del spline natural evaluadas en los valores de los datos. El vector y contiene los valores de respuesta y1,…,yn. En la práctica, existen algoritmos muy eficientes para calcular \hat {\beta }\) en (4) . En lugar de especificar una base spline natural para f, también es posible trabajar con una base B-spline sin restricciones, ya que la penalización en (3) impone automáticamente las restricciones de linealidad en los nodos x(1) y x(n) (ver , Capítulo 5, y , Capítulo 2). En lo que respecta a la base B-spline, los resultados de la estimación no dependerán de la elección de los nodos límite: es posible utilizar x(1) y x(n) como nodos límite o incluir x(1) y x(n) en el conjunto de nodos interiores.
Si n es grande y el intervalo está cubierto densamente por los datos observados, normalmente no es necesario colocar un nudo en cada xi,i=1,…,n. En su lugar, el spline de suavizado puede ser aproximado por un spline de regresión penalizado que utiliza un conjunto reducido de nodos. Una clase muy popular de splines de regresión penalizados son los P-splines , que se basan en la base cúbica B-spline y en un conjunto «grande» de nodos equidistantes (normalmente, 10-40). En lugar de evaluar la integral en (3), las P-splines se basan en una penalización de segundo orden definida por
que, en caso de nudos uniformemente espaciados, puede demostrarse que es una aproximación a Jβ. El operador diferencial de segundo orden Δ2 se define por Δ2βk:=(βk-βk-1)-(βk-1-βk-2). Por tanto, la penalización puede expresarse como β⊤Pβ, donde P se define por D⊤D con D una matriz de diferencias. Se deriva fácilmente que el estimador resultante de β tiene la misma estructura que 2, con Ω sustituido por P.
Una propiedad conveniente de las P-splines es que son numéricamente estables y muy fáciles de definir e implementar. En particular, es mucho más fácil configurar la matriz de diferencias D que la matriz Ω. Además, es sencillo extender la penalización Jβ (y, por tanto, la matriz D) a diferencias de orden superior Δq con q>2. También es posible utilizar una secuencia de nudos que no esté uniformemente espaciada; en este caso, es necesario introducir pesos. Debido a que las P-splines con nudos desigualmente espaciados rara vez se utilizan en la práctica, no las consideramos aquí y nos referimos a ellas en su lugar.
Las splines de suavizado y las P-splines superan el problema de la selección de nudos hasta cierto punto. Su filosofía es utilizar un gran número de nudos y luego dejar que λ controle la cantidad de suavidad. Esto da lugar a un parámetro de ajuste adicional, sin que exista un consenso general sobre cómo ajustar este parámetro. Algunas formas populares de determinar el valor «óptimo» de λ utilizan la validación cruzada generalizada (GCV), el AIC o una representación de modelo mixto.
Splines en R
El paquete de instalación básico de R contiene un conjunto de funciones que pueden ajustar splines polinomiales simples y splines de suavizado. Se incluyen más funciones en la biblioteca splines escrita por DM Bates y WN Venables. El paquete ha sido el caballo de batalla del ajuste de splines durante muchos años y ahora forma parte de la distribución básica de R. Hay más de 100 otros paquetes que dependen de splines al cargarse. El paquete contiene varias funciones para crear bases de splines, como bs para B-splines y ns para splines naturales, que son ampliamente utilizadas, pero también algunas funciones más especializadas para crear funciones base (como periodicSpline que crea una interpolación periódica de splines) o comandos que son útiles como el comando predict.bSpline que evaluaría una spline en los nuevos valores de X.
Los valores por defecto de bs crearán una base B-spline cúbica con dos nudos de frontera y uno interior situados en la mediana de los valores de los datos observados. El usuario puede conseguir una mayor flexibilidad, aumentando la colocación y el número de nudos y/o cambiando su ubicación. La figura 3 (código 2 en el archivo adicional 1: Apéndice) muestra las B-splines creadas con diferentes opciones. La parte superior presenta splines lineales, es decir, polinomios de primer orden (el grado es uno) conectados entre sí en nudos equidistantes. La parte inferior presenta polinomios cúbicos (grado 3).
Hay que tener en cuenta que los B-splines creados en R con bs() están automáticamente acotados por el rango de los datos, y que los nodos adicionales (τ1,…,τd) se establecen igual a los nodos de frontera, dando múltiples nodos en ambos extremos del dominio. Este enfoque es útil en casos univariantes y tiene algunas características atractivas desde el punto de vista computacional. Sin embargo, si se trabaja en un problema de suavizado bidimensional, utilizando productos tensoriales de B-splines, o cuando se trabaja con P-splines, esta base es inadecuada y puede conducir a resultados espurios.
Los splines naturales se pueden crear dentro del paquete splines, utilizando el comando ns. Por defecto, a menos que el usuario especifique los grados de libertad o los nodos, la función devuelve una línea recta dentro de los nodos límite. La figura 4 (código 3 en el archivo adicional 1: Apéndice muestra splines naturales creados con diferentes opciones.
Para ilustrar cómo se pueden utilizar estas funciones en la práctica, considere de nuevo los datos de la sección 2.0.1. La figura 5 (creada por (código 4 en el archivo adicional 1: Apéndice)) muestra los ajustes obtenidos utilizando los siguientes comandos: poly() para splines polinómicos ortogonales simples, smooth.spline() para splines de suavizado, bs() y ns() de la biblioteca splines, para splines B y splines naturales respectivamente. El gráfico superior izquierdo muestra un ajuste lineal simple sobre los datos (línea discontinua) y un ajuste polinómico de tercer grado que es capaz de capturar la relación más compleja entre las variables. Sin embargo, el gráfico de la esquina superior derecha es especialmente interesante, ya que presenta los ajustes utilizando los valores por defecto de las funciones spline. La línea verde proviene de las funciones poly() y ns() que por defecto definen una línea recta. En el otro extremo, la línea azul es un ajuste de la función smooth.spline() que, si no se especifican grados de libertad, tiende a suavizar los datos, es decir, produce un ajuste muy flexible basado -aquí- en 45 grados de libertad. Se puede conseguir un ajuste -visualmente- razonable de los datos cuando se especifican cuatro grados de libertad (gráfico inferior izquierdo). Se puede ver que hay algunas diferencias dependiendo de la base elegida. La base polinómica (línea negra) es un poco más flexible que el resto, especialmente a edades más altas. Por otro lado, un spline de suavizado restringido a sólo cuatro grados de libertad es más rígido que los otros enfoques, pero probablemente suaviza en exceso los datos a edades pequeñas, entre los años 0 y 10. Entre los dos extremos, los splines B y los splines naturales proporcionan ajustes muy similares que capturan el efecto de las edades pequeñas y tienden a estar menos influenciados por los casos extremos al final del espectro de edad. Por último, el gráfico inferior derecho muestra cómo los ajustes se vuelven mucho más flexibles con grados de libertad adicionales y sugiere un potencial sesgo de sobreajuste debido al uso de grados de libertad excesivos.
Una nota sobre los grados de libertad
En la práctica, siempre es útil definir un spline por grados de libertad. Este enfoque es particularmente útil cuando se trabaja con B-splines y splines naturales. Los B-splines tienen d+K, mientras que una función base de un spline cúbico natural con K nudos tiene K+1 grados de libertad, respectivamente. Por defecto, la función bs en R crea B-splines de grado 3 sin nodos interiores y nodos de frontera definidos en el rango de la variable X. Como tal, la función crea tres funciones de base. Consideremos ahora el siguiente caso: cuando un usuario define una B-spline con un nudo interior en la mediana de X (bs(x,nudos=mediana(x))) el software creará cuatro funciones (d=3 más K=1 nudos interiores, cuatro grados de libertad). Sin embargo, si el usuario especifica en la función los nudos interiores dentro del argumento nudos (bs(x,nudos=c(min(x),mediana(x),max(x)))), la función tendrá seis grados de libertad (d=3 más k=3). Se debe tener una precaución similar con la función ns.
Cuando se trabaja con splines de suavizado, no es fácil especificar los grados de libertad, ya que variarán dependiendo del tamaño de la penalización. Sin embargo, en la práctica, los splines penalizados también pueden restringirse a un número máximo de grados de libertad o a los grados de libertad deseados.
Otros paquetes de splines
En términos generales, la lista ampliada de paquetes de splines contiene o bien enfoques bastante similares a lo que se presenta aquí o bien casos muy especializados que se dirigen a aplicaciones específicas. En la Tabla 1 se presentan algunos de estos paquetes junto con el número de descargas. El número se refiere al número de veces que se ha descargado un paquete, pero no a los usuarios únicos. Está fuera del alcance de este trabajo describir en detalle todos estos enfoques.