Definición.

Cuando se observa que una variable presenta diferente comportamiento en diversos tramos de su dominio, un método usual es dividirlo en segmentos y ajustar la curva adecuada en cada segmento. Las funciones spline ofrecen una forma sencilla, flexible y útil para realizar un ajuste por segmentos.

Los son polinomios de orden \(k\) por segmentos. Los puntos de unión de los segmentos se llaman “nudos” (knots) generalmente se requiere que los valores de la función y de las primeras \(k-1\) derivadas concuerden con los nudos, para que el sea una función continua con \(k-1\) derivadas continuas. La cúbica suele ser muy adecuada para muchos problemas prácticos.

Una función cúbica con \(h\) nudos, \(t_{1}<t_{2}<t_{3}<\cdots t_{h}\), con primeras y segundas derivadas continuas se puede expresar: \[\begin{equation*} E(y)=S(x)=\displaystyle\sum_{j=0}^{3}\beta_{0}x^{j}+\displaystyle\sum_{i=1}^{h}\beta_{i}(x-t_{i})^{*} \end{equation*}\] en donde:

\[\begin{equation*} (x-t_{i})^{*}= \begin{cases} (x-t_{i}) &\text{si $x-t_{i}>0$}\\ 0 & \text{si $x-t_{i}\leq 0$} \end{cases} \end{equation*}\]

Se parte de que se conocen las posiciones de los nudos, pero si ellos son parámetros que se deben estimar, el problema resultante lleva a una regresión no lineal; sin embargo, si se conocen las posiciones de los nudos, la ecuación puede ser estimada aplicando mínimos cuadrados ordinarios en forma directa.

La decisión que debe tomarse a priori es la ubicación de los nudos y el grado del polinomio por ajustar. Sin embargo, algunos sugieren que debe tenerse un mínimo posible de nudos y en cada uno de los tramos definidos contar con al menos cinco puntos muestrales, y, hasta donde sea posible, los puntos extremos deberían estar centrados en el segmento y los puntos de inflexión deberían estar cerca de los nudos.

El modelo de spline cúbico se puede expresar de la forma:

\[\begin{equation*} y=\beta_{0}+\beta_{1}x+\beta_{2}x^{2}+\beta_{3}x^{3}+\beta_{4}(x-n_{1})^{3}+\cdots \beta_{n}(x-n_{i})^{3} \end{equation*}\]

donde

\(n_{1}\cdots n_{i}\) son los \(i\) nudos

Ejemplo

Se utiliza el ejemplo del libro “Introducción al análisis de regresión lineal” de Montgomery-Peck-Vining, página 206.

Para el modelo se utilizan como nudo los valores 6.5 y 13.

x<-seq(0,20,0.5)
y<-c(8.33,8.23,7.17,7.14,7.31,7.60,7.94,8.30,8.76,8.71,9.71,10.26,10.91,11.67,
11.76,12.81,13.30,13.88,14.59,14.05,14.48,14.92,14.37,14.63,15.18,14.51,14.34,
13.81,13.79,13.05,13.04,12.60,12.05,11.15,11.15,10.14,10.08,9.78,9.80,9.95,9.51)
x2<-x^2
x3<-x^3
z<-(x-6.5)^3
x4<-ifelse(z<0,0, z)
w<-(x-13)^3
x5<-ifelse(w<0,0,w)
B<-data.frame(x,x2,x3,x4,x5,y)
head(B)
##     x   x2     x3 x4 x5    y
## 1 0.0 0.00  0.000  0  0 8.33
## 2 0.5 0.25  0.125  0  0 8.23
## 3 1.0 1.00  1.000  0  0 7.17
## 4 1.5 2.25  3.375  0  0 7.14
## 5 2.0 4.00  8.000  0  0 7.31
## 6 2.5 6.25 15.625  0  0 7.60
modelo<-lm(y~x+x2+x3+x4+x5)
summary(modelo)
## 
## Call:
## lm(formula = y ~ x + x2 + x3 + x4 + x5)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.45168 -0.18499 -0.03547  0.20577  0.61694 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8.465678   0.200520  42.219  < 2e-16 ***
## x           -1.453124   0.181586  -8.002 2.04e-09 ***
## x2           0.489889   0.043018  11.388 2.54e-13 ***
## x3          -0.029467   0.002848 -10.347 3.44e-12 ***
## x4           0.024706   0.004039   6.116 5.43e-07 ***
## x5           0.027112   0.003578   7.577 6.98e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2678 on 35 degrees of freedom
## Multiple R-squared:  0.9904, Adjusted R-squared:  0.9891 
## F-statistic: 725.5 on 5 and 35 DF,  p-value: < 2.2e-16

\[\begin{equation*} y= 8.465678-1.453124x+0.489889x^{2}-0.029467x^{3}+0.024706(x-6.5)^{3}+\\0.027112(x-13)^{3} \end{equation*}\]

Se obtienen los valores estimados por el modelo:

pre1<-predict(modelo)
pre1
##         1         2         3         4         5         6         7         8 
##  8.465678  7.857905  7.472976  7.288791  7.283249  7.434250  7.719694  8.117480 
##         9        10        11        12        13        14        15        16 
##  8.605508  9.161678  9.763889 10.390041 11.018034 11.625768 12.194229 12.716760 
##        17        18        19        20        21        22        23        24 
## 13.189788 13.609745 13.973058 14.276156 14.515470 14.687427 14.788458 14.814991 
##        25        26        27        28        29        30        31        32 
## 14.763456 14.630282 14.411897 14.108121 13.732326 13.301276 12.831735 12.340464 
##        33        34        35        36        37        38        39        40 
## 11.844228 11.359789 10.903909 10.493353 10.144883  9.875262  9.701253  9.639619 
##        41 
##  9.707124

El gráfico con los datos originales y los valores estimados se muestran a continuación.

plot(x,y,pch=16)
lines(x,pre1,lwd=2,col="red")

Utilizando el paquete splines

library(splines)
modelo2<-lm(y~bs(x ,knots =c(6.5,13) ))
summary(modelo2)
## 
## Call:
## lm(formula = y ~ bs(x, knots = c(6.5, 13)))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.45168 -0.18499 -0.03547  0.20577  0.61694 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                  8.4657     0.2005  42.219  < 2e-16 ***
## bs(x, knots = c(6.5, 13))1  -3.1484     0.3934  -8.002 2.04e-09 ***
## bs(x, knots = c(6.5, 13))2   4.3532     0.2843  15.312  < 2e-16 ***
## bs(x, knots = c(6.5, 13))3   8.5518     0.3691  23.169  < 2e-16 ***
## bs(x, knots = c(6.5, 13))4   0.5990     0.3059   1.958 0.058192 .  
## bs(x, knots = c(6.5, 13))5   1.2414     0.2871   4.324 0.000121 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2678 on 35 degrees of freedom
## Multiple R-squared:  0.9904, Adjusted R-squared:  0.9891 
## F-statistic: 725.5 on 5 and 35 DF,  p-value: < 2.2e-16

Los valores predichos por este modelo

pre2<-predict(modelo2)
pre2
##         1         2         3         4         5         6         7         8 
##  8.465678  7.857905  7.472976  7.288791  7.283249  7.434250  7.719694  8.117480 
##         9        10        11        12        13        14        15        16 
##  8.605508  9.161678  9.763889 10.390041 11.018034 11.625768 12.194229 12.716760 
##        17        18        19        20        21        22        23        24 
## 13.189788 13.609745 13.973058 14.276156 14.515470 14.687427 14.788458 14.814991 
##        25        26        27        28        29        30        31        32 
## 14.763456 14.630282 14.411897 14.108121 13.732326 13.301276 12.831735 12.340464 
##        33        34        35        36        37        38        39        40 
## 11.844228 11.359789 10.903909 10.493353 10.144883  9.875262  9.701253  9.639619 
##        41 
##  9.707124

gráfico del modelo

plot(x,y,pch=16)
lines(x,pre2,lwd=2,col="red")

Error cuadrático medio (ecm) y error medio absoluto (ema) de los dos procedimientos

ecm1<- mean((y - pre1)^2);ecm1
## [1] 0.06122573
ema1<-mean(abs(y-pre1));ema1
## [1] 0.2093843
ecm2<- mean((y - pre2)^2);ecm2
## [1] 0.06122573
ema2<-mean(abs(y-pre2));ema2
## [1] 0.2093843

.