Regression splines

Questo metodo per fare regressione è non parametrico e si è costruito per risolvere il problema dell’elevata collinearità presente in molti modelli polinomiali; infatti se ad esempio supponiamo:

\[y = f(x) + \varepsilon\]

Consideriamo la seguente funzione per \(x\) \[f(x) = sin(2(4x-2)) + 2e^{-(16^2)(x-.5)^2}\]

Servirebbe un elevato grado (pari a 15) alla regressione polinomiale per descrivere in modo accurato la funzione ma questo comporterebbe un alta correlazione che ha come consueguenze:

• Instabilità delle stime
• Sensibilità agli outlier

L’idea che sta alla base di questo modello è quello di dividere lo spazio da analizzare in più sezioni come divisione si scelgono dei punti detti nodi (knots) attraverso un qualche criterio (solitamente i quantili) e in ogni sezione si attua una regressione.
Utilizzando questo metodo come sopra descritto otteniamo delle funzioni a tratti (piecewise):

• Definiamo \(K\) knots (internal breakpoints) nel range di valori di \(x\): \[\min(x)<\xi_1< \ldots < \xi_K<\max(x)\]

• poi dobbiamo fittare il modello polinomiale con d gradi nei \(K+1\) intervalli

\[(-\infty, \xi_{1}], (\xi_{1}, \xi_{2}], \ldots, (\xi_{K-1}, \xi_{K}], (\xi_{K}, +\infty)\]

Consideriamo quindi una funzione generale per \(x\):

\[f(x) = \sum_{j=1}^{p} \beta_j b_j(x)\] dove \[b_j(\cdot)\] sono funzioni note chiamate basis functions.

set.seed(123)
x = sort(runif(500))
mu = sin(2*(4*x-2)) + 2*exp(-(16^2)*((x-.5)^2))
y = rnorm(500, mu, .3)
knots = seq(0.1, 0.9, by=.1)
xcut = cut(x, c(-Inf,knots,Inf) )

Step function d=0 \[b_1(X) = I\{x < \xi_1\}, b_2(x) = I\{\xi_1 \leq x < \xi_2\} ,\ldots, b_{K}(x) = I\{\xi_{K-1} \leq x < \xi_K\}, b_{K+1}(x) = I\{ x \geq \xi_K\}\]

plot(x,y ,col="lightgray")
lines(x, mu, col=2, lwd=2)
abline(v=knots, lty=2)
for (i in 1:length(levels(xcut)) ){
sub = (xcut==levels(xcut)[i])
lines(x[sub], fitted(lm(y[sub] ~ 1)) )
}

Piecewise linear regression d=1

\[b_1(x) =I\{x < \xi_1\}, b_2(x) = x \cdot I\{x < \xi_1\}, b_3(x) = I\{\xi_1 \leq x < \xi_2\}, b_4(x) = x \cdot I\{\xi_1 \leq x < \xi_2\},\ldots, b_{2(K+1)}(x) = x \cdot I\{ x \geq \xi_K\}\]

Piecewise polynomial (cubic) regression d=3 \[b_1(x) = 1, b_2(x)= x, b_{3}(x) = x^2, b_{4}(x) = x^3\]

plot(x,y ,col="lightgray")
lines(x, mu, col=2, lwd=2)
abline(v=knots, lty=2)
for (i in 1:length(levels(xcut)) ){
sub = (xcut==levels(xcut)[i])
lines(x[sub], fitted(lm(y[sub] ~ poly(x[sub],3))) )
}

In questa prima soluzione, come possiamo notare, però non è garantita la continuità. Per risolvere questo problema si è pensato di imporre che in ogni nodo la funzione abbia derivate continue dal grado 0 al grado d-1:
VINCOLI DI CONTINUITA’
\[f(\xi_i) = Y_i \ \ \ \ \ \ \ \ con \ \ \ \ \ \ \ i=1,...,K \] \[f(\xi_i^{+}) = f(\xi_i^{-}) \ \ \ \ \ \ \ \ con \ \ \ \ \ \ \ i=2,...,Z-1\] \[f'(\xi_i^{+}) = f'(\xi_i^{-}) \ \ \ \ \ \ \ \ con \ \ \ \ \ \ \ i=2,...,K-1\] \[f''(\xi_i^{+}) = f''(\xi_i^{-}) \ \ \ \ \ \ \ \ con \ \ \ \ \ \ \ i=2,...,K-1\] \[.\]\[.\]\[.\]

Trade off tra gradi e nodi

Qui sotto c'è una applicazione che mostra come i due paratri influenzano lo sviluppo dell'algoritmo (questa applicazione è stata fatta con la versione graduita di shinyapp quindi mi scuso per eventuali disagi fornisco anche il codice .Rmd cosi da poter agevolare la visualizzazione: codice.Rmd)

Questo ragionamento per quanto corretto pone due problemi:

• Oltre la derivata seconda risulta difficile valutare la continuità di una funzione.
• I parametri da stimare sono \((d+1)(K-1)\) mentre i vincoli sono soltanto \(K+(d(K-d+1))\)

Per quanto riguarda il primo problema si utilizza quasi universalmente un grado massimo pari a \(d=3\) e si parla di Natural cubic splines; invece per risolvere il fatto che non si identifica una univoca funzione solitamente si impongono che le derivate seconde negli intervalli estremi siano nulle.

In definitiva la natural cubic spline con \(K\) knots \(\xi_1,\ldots,\xi_K\) è una funzione lineare a tratti \(f\) che:
\(f\) è un polinomio di grado 3 \[[\xi_1,\xi_2],[\xi_2,\xi_3],\ldots,[\xi_{K-1},\xi_{K}]\]
\(f\) è lineare negli intervalli \((-\infty, \xi_1]\) e \([\xi_K,\infty)\)
\(f\), \(f'\) e \(f''\) sono continue in ogni nodo \(\xi_1,\ldots,\xi_K\): \(f(\xi_k^-)=f(\xi_k^+)\), \(f'(\xi_k^-)=f'(\xi_k^+)\) e \(f''(\xi_k^-)=f''(\xi_k^+)\), \(k=1,\ldots,K\), dove \(\xi_k^+\) e \(\xi_k^-\) indicano il limite destro e sinistro della funzione \(f(\cdot)\) nel nodo \(\xi_k\)

Numero di nodi

Possiamo evitare del tutto il problema della selezione del nodo formulando un problema di minimizzazione penalizzato:

\[\underset{f \in \mathcal{F}''} \min \,\, \left\{ \sum_{i=1}^{n}(y_i - f(x_i))^2 + \lambda \int \{f''(t)\}^2 dt \right\}\] dove \(f\) si distribuisce come \(\mathcal{F}''\) of funzioni due volte differenziabili

• Il primo termine è RSS, e confronta \(f(x_i)\) con \(y_i\) in ogni \(x_i\)

• Il secondo termine è una roughness penalty (rugosità dela funzione) e controllo come fluttua (wiggly) come si muove \(f(x)\)

• La penalità dipende dal tuning parameter \(\lambda\geq 0\):
• \(\lambda=0\) non impone condizione su \(f\) e viene calcolata interpolando con i dati

• \(\lambda = \infty\) rende la curvatura impossibili, e quindi la funzione \(f\) sarà lineare

Teorema
Per risolvere questo problema di minimizzazione penalizzata \(f\) tra tutte le funzioni possibili di \(\mathcal{F}''\) deve essere una natural cubic spline

• Teorema (Green and Silverman, 1994)

Out of all twice-differentiable functions \(f\), the one that minimizes \[\sum_{i=1}^{n}(y_i - f(x_i))^2 + \lambda \int \{f''(t)\}^2 dt\] is a natural cubic spline with knots at every unique value of \(x_i\)