Il termine ‘spline’ si riferisce a uno strumento artigianale, una sottile striscia flessibile di legno o metallo, usata per disegnare curve lisce. Diversi pesi verrebbero applicati su varie posizioni in modo che la striscia si pieghi secondo il loro numero e posizione. Questa sarebbe costretta a passare attraverso una serie di punti fissi: perni metallici, le costole di una barca, ecc. Su una superficie piana questi erano spesso dei pesi con un gancio attaccato e quindi facili da manipolare. La forma del materiale piegato prenderebbe naturalmente la forma di una curva spline. Allo stesso modo, le spline sono utilizzate in statistica per riprodurre matematicamente forme flessibili. I nodi sono posti in diversi punti all’interno della gamma di dati, per identificare i punti in cui pezzi funzionali adiacenti si uniscono tra loro. Invece di strisce di metallo o legno, si scelgono pezzi funzionali lisci (di solito polinomi di basso ordine) per adattare i dati tra due nodi consecutivi. Il tipo di polinomio e il numero e il posizionamento dei nodi è ciò che definisce il tipo di spline.
Esempio motivante
Con l’introduzione dei modelli additivi generalizzati (GAM) nel 1986, l’uso della modellazione spline è diventato uno strumento consolidato nell’analisi di regressione statistica. Per illustrare ciò, consideriamo i dati su un set di 892 femmine sotto i 50 anni raccolti in tre villaggi dell’Africa occidentale (dati disponibili nel file aggiuntivo 1: Appendice). Vorremmo esplorare la relazione tra l’età (in anni) e una misura grezza del grasso corporeo, che è lo spessore della plica tricipitale. La figura 1 mostra la relazione tra l’età e lo spessore della plica tricipitale misurata in scala logaritmica. Per ulteriori informazioni sui dati vedere .
Un semplice modello di regressione della forma yi=β0+β1xi+ε,i=1,…,n, difficilmente darebbe un’approssimazione del modello osservato, poiché è ovvio che la relazione non è lineare. Il modello può essere esteso per accogliere effetti non lineari usando alcuni polinomi. Quindi, gli effetti non lineari potrebbero essere modellati da un polinomio di grado 3 dato da:
dove u è una funzione di x chiamata funzione base, definita qui da:
Il modello di regressione descritto nella Eq. 1 è ancora un modello lineare, nonostante il fatto che fornisce una funzione non lineare della variabile predittiva. Il modello è ancora lineare nei coefficienti e può essere adattato usando i metodi ordinari dei minimi quadrati. La base può essere creata in R usando la funzione poly(x,3) con gli input x (riferito alla variabile), e p (riferito al grado del polinomio). Questo porta a un semplice modello univariato liscio della forma: yi=f(xi)+ε dove f() è qualche funzione/trasformazione del predittore. Tale modello può essere facilmente adattato in R usando: lm(y ∼poly(x,3)). Nonostante la semplicità, la regressione polinomiale ha diversi svantaggi, il più importante dei quali è la non località. Ciò significa che la funzione adattata a un dato valore x0 dipende dai valori dei dati lontani da quel punto. È facile vedere questo in azione adattando una polinomiale a un insieme di dati e spostando uno dei punti dati vicino al bordo destro verso l’alto o verso il basso. Come risultato, la funzione adattata di solito cambierà lontano da quella coordinata x.
Considera, invece di adattare un polinomio globale, di partizionare l’intervallo di x in intervalli più piccoli, utilizzando un numero arbitrario e una posizione di punti, τ, chiamati anche nodi. Un semplice modello continuo piecewise può essere montato definendo le funzioni: f1(x)=1,f2(x)=x,f3(x)=(x-τ1)+,f4(x)=(x-τ2)+,…, con “+” una funzione definita come:
L’insieme di queste funzioni porta ad una funzione composita f(x).
Definizione delle spline
La spline metallica del disegnatore può assumere forme arbitrarie, per esempio la sezione trasversale dell’ala di un aereo o la spirale di una pompa centrifuga. Per le applicazioni statistiche assumeremo curve della forma f(X), cioè un singolo valore y per ogni x. Il predittore x può essere una singola variabile o variabili multiple. La nostra discussione si concentrerà quasi interamente su una funzione univariata con \(X\in \mathbb {R}\). Definiamo un insieme di nodi τ1<…<τK nell’intervallo di X. Una spline f(X) sarà una funzione liscia, che soddisfa certe proprietà di differenziabilità menzionate più avanti, tale che f(X) è un polinomio di grado d. Le spline di legno o metallo hanno derivate continue di tutti gli ordini, poiché sono un oggetto fisico. Questo non è vero per le spline statistiche. Piuttosto imponiamo un criterio di morbidezza per cui tutte le derivate di ordine inferiore a d sono continue. Una spline fisica è lineare oltre l’ultimo nodo e possiamo imporre un ulteriore vincolo: le derivate di ordine 2 o maggiore sono nulle nei nodi più a sinistra e più a destra; le spline con questo ulteriore vincolo sono note come spline “ristrette” o “naturali”. Per ottenere curve più flessibili si può aumentare il numero di nodi o il grado del polinomio. C’è comunque un compromesso: aumentare il numero di nodi può sovraadattare i dati e aumentare la varianza, mentre diminuire il numero di nodi può risultare in una funzione rigida e restrittiva che ha più bias.
Rappresentazione tramite funzioni base
Assumiamo che la funzione incognita f sia rappresentata da una funzione spline con sequenza di nodi fissa e grado d fisso. Poiché queste ultime funzioni formano uno spazio vettoriale V, è possibile scrivere f come
dove le Bk sono un insieme di funzioni base che definiscono V e βk sono i coefficienti spline associati. Con k nodi ci sono k+1 polinomi di grado d insieme a d∗k vincoli, portando a (d+1)(k+1)-d∗k=d+k+1 parametri liberi; per una spline naturale ci sono k parametri liberi. Poiché βB=(βA)(A-1B)=γB∗ per qualsiasi matrice non singolare A ci sono un numero infinito di possibili set di basi per l’adattamento della spline.
La rappresentazione in (2) ha il vantaggio che la stima di f si riduce alla stima dei coefficienti βk. Più specificamente, l’espressione in (2) è lineare nel vettore dei coefficienti β=(β1,…,βK+d+1). Quindi la stima di f può essere vista come un problema di ottimizzazione che è lineare nelle variabili trasformate B1(X),…,BK+d+1(X), permettendo l’uso di tecniche di stima ben consolidate per l’uso di spline in una vasta gamma di modelli di regressione multivariabile (generalizzata). È importante notare che la modellazione delle spline riduce la stima delle funzioni f() alla stima di un piccolo insieme di coefficienti con valore reale.
Come sottolineato da vari autori (ad esempio, l’alta flessibilità della modellazione delle spline viene al prezzo di un certo numero di parametri di regolazione. Due di questi, la scelta delle funzioni di base B e il grado d dei polinomi sottostanti risultano avere un impatto minimo. Infatti, gli adattamenti delle spline sono notevolmente robusti al grado d. I polinomi cubici (d=3) sono lo standard abituale in quanto risultano in curve che appaiono perfettamente lisce all’occhio umano. Se le derivate delle curve adattate sono di interesse, un ordine superiore è talvolta appropriato, ma in generale gli adattamenti per d>3 sono effettivamente indistinguibili. I fits con d=1 o d=2 hanno proprietà statistiche quasi identiche, ma appaiono più frastagliati. La scelta tra due set di basi B e B∗ per definizione non cambierà le previsioni di un fit e quindi si riduce a questioni di convenienza.
Le due scelte chiave sono il numero e la spaziatura dei nodi e l’uso (o meno) di una funzione di penalità, ad esempio, la derivata seconda integrata della spline. Quando non c’è penalità, la creazione delle variabili trasformate può essere fatta separatamente e le nuove variabili sono semplicemente incluse in un modello standard; non è richiesta alcuna modifica della procedura di regressione sottostante. Questo approccio è spesso indicato come spline di regressione; la flessibilità della funzione non lineare risultante è interamente una funzione del numero di nodi. L’inclusione di una penalità di smoothing, d’altra parte, richiede la modifica della routine di adattamento per poterla ospitare. Questa deve essere inclusa in ogni funzione di regressione separatamente. Le spline di lisciatura risultanti hanno diverse proprietà desiderabili, ma la complessità aggiunta della funzione di lisciatura può essere una ragione per non essere stata usata più spesso in contesti applicativi.
Anche se sono state condotte notevoli ricerche per esplorare le proprietà matematiche dei vari approcci spline (vedi , gli statistici applicati e gli analisti di dati difficilmente sembrano essere consapevoli di questi risultati quando usano la modellazione spline in applicazioni pratiche. Infatti, molti degli articoli identificati dalla nostra ricerca sul web non contengono alcuna giustificazione sulla logica per la scelta del metodo spline utilizzato.
Basi spline popolari
Ci sono numerose opzioni per la definizione delle funzioni di base Bk, dove le varie basi spline differiscono rispetto alle loro proprietà numeriche. In questa sezione, introdurremo alcune delle basi spline più popolari, vale a dire la base della serie di potenza troncata, la base B-spline e la base spline cardinale.
Serie di potenza troncata e spline cubiche
La base della serie di potenza troncata è definita dalle funzioni base
Un vantaggio delle funzioni base di cui sopra è la loro facile interpretazione: Partendo da un polinomio “base” di grado d definito su (prima linea di equazione), le deviazioni dal polinomio base sono successivamente aggiunte alla funzione spline a destra di ciascuno dei K nodi (seconda linea). Una spline troncata a base di potenza è d-1 volte differenziabile ai nodi e ha d+K gradi di libertà. È relativamente facile per l’utente creare una serie di potenza troncata in R. Lasciamo che x rappresenti alcune osservazioni in , allora una base di potenza troncata di grado d=3 con 5 nodi equamente distanziati lungo l’intervallo di x può essere creata utilizzando il codice 1 nel file aggiuntivo 1: Appendice (Fig. 2).
Una caratteristica delle serie di potenza troncate è che i supporti delle funzioni non sono locali, con alcuni dei Bk definiti sull’intero intervallo di dati. Questo potrebbe portare ad alte correlazioni tra alcune spline di base, implicando instabilità numeriche nella stima delle spline. Per la base della serie di potenza troncata, un esempio è dato in , Capitolo 5.
Le spline cubiche sono create utilizzando un polinomio cubico in un intervallo tra due nodi successivi. La spline ha quattro parametri su ciascuna delle regioni K+1 meno tre vincoli per ogni nodo, risultando in un K+4 gradi di libertà.
Una funzione spline cubica, con tre nodi (τ1,τ2,τ3) avrà 7 gradi di libertà. Usando la rappresentazione data nell’Eq. 2, la funzione può essere scritta come:
B-spline
La base B-spline è una base spline comunemente usata che è basata su una speciale parametrizzazione di una spline cubica. La base B-spline, è basata sulla sequenza di nodi
dove gli insiemi ξd+2 := τ1,…,ξd+K+1:=τK e ξd+1:=a,ξd+K+2:=b sono chiamati rispettivamente “nodi interni” e “nodi di confine”. La scelta dei nodi aggiuntivi ξ1,…,ξd e ξd+K+3,…,ξ2d+K+2 è essenzialmente arbitraria. Una strategia comune è quella di porli uguali ai nodi di confine. In alternativa, se i nodi interni e i nodi di confine ξd+1<…<ξd+K+2 sono scelti per essere equidistanti, cioè, ξk+1-ξk=δ ∀k∈{d+1,…,d+K+1}, i nodi limite possono essere posti a ξd+1-δ,…,ξd+1-d-δ e ξd+K+2+δ,…,ξd+K+2+d-δ.
Per d>0, le funzioni base B-spline di grado d (indicate con \(B_{k}^{d}(x)\)) sono definite dalla formula ricorsivaFootnote 1
dove
e \(B_{k}^{0}(x) \equiv 0\) se ξk=ξk+1. Le B-spline hanno il vantaggio che le funzioni base hanno un supporto locale. Più specificamente, sono più grandi di zero negli intervalli attraversati da d+2 nodi e zero altrove. Questa proprietà risulta in un’alta stabilità numerica, e anche in un algoritmo efficiente per la costruzione delle funzioni di base, vedi per i dettagli.
Splines naturali cubiche e cardinali
Una spline polinomiale come una cubica o una B-spline, può essere erratica ai confini dei dati. Per affrontare questo problema, le spline naturali sono spline cubiche che hanno il vincolo aggiuntivo di essere lineari nelle code dei nodi di confine (-∞,a],, Capitolo 4.
Oltre alle spline naturali in serie di potenza troncata, B-spline e spline cardinali, esistono varie altre basi – meno popolari. Per una panoramica, facciamo riferimento ai libri di .
Spline penalizzate
Le spline presentate finora sono spesso chiamate spline di regressione. Oltre alla scelta della base della spline (B-spline, serie di potenza troncata, ecc.), bisogna scegliere il numero di nodi e le posizioni dei nodi. Ovviamente, questi parametri di regolazione possono avere un impatto importante sulla forma stimata di una funzione spline: Un grande numero di nodi implica un’alta flessibilità, ma può anche risultare in un overfitting dei dati a disposizione. Al contrario, un piccolo numero di nodi può risultare in una stima “oversmooth” che è incline a bias di under-fit (vedi ).
Un approccio popolare per facilitare la scelta delle posizioni dei nodi nella modellazione spline è l’uso di spline penalizzate. Dato un campione i.i.d. di dati (x1,y1),…(xn,yn), una spline penalizzata è la soluzione del problema
dove lβ denota la log-likelihood (o, nel caso della regressione di Cox, la log-liquidità parziale) e Jr è una penalità di rugosità che diventa piccola se la funzione spline è “liscia”. Generalmente, le spline penalizzate si basano sull’idea che la funzione incognita f sia modellata da una spline con un gran numero di nodi, permettendo un alto grado di flessibilità. D’altra parte, una stima approssimativa della spline che ha un alto valore di lβ ed è vicina ai valori dei dati risulta in un grande valore di Jβ. La massimizzazione di questa funzione implica quindi un trade-off tra la morbidezza e l’adattamento del modello che è controllato dal parametro di sintonizzazione λ≥0.
Un caso speciale è il problema dei minimi quadrati penalizzati
in regressione gaussiana. La penalità \(J_{\beta } \, \int _{a}^{b} \sinistra (\parziale ^{2} f / \parziale x^{2} destra)^{2} dx\) esprime la “scorrevolezza” di una funzione spline in termini di derivata seconda di f. Per data λ, si può dimostrare che la soluzione è una spline cubica naturale con sequenza di nodi x(1)<…<x(n), cioè, le posizioni dei nodi non devono essere scelte ma sono “naturalmente” date dai valori unici ordinati dei dati di X. In letteratura, questo tipo di spline è indicato come spline di smoothing. Da notare che si può dimostrare che una spline di lisciatura interpola i dati se λ=0, mentre λ=∞ implica una funzione lineare. Si noti che le spline di lisciatura sono un caso speciale della classe più generale delle spline a piastra sottile, che permettono un’estensione del criterio in Eq. (3) a xi più grandi (vedi , Sezione 4.15], e per i dettagli).
Una proprietà conveniente delle spline di lisciatura è che la penalità Jβ può essere scritta come β⊤Ωβ con una matrice di penalità opportunamente definita Ω. Pertanto la soluzione della (3) è data dalla stima dei minimi quadrati penalizzati
dove B è una matrice di dimensione n×n contenente le funzioni base spline naturali valutate ai valori dei dati. Il vettore y contiene i valori di risposta y1,…,yn. In pratica, esistono algoritmi molto efficienti per calcolare \(\che {\beta}) in (4) . Invece di specificare una base spline naturale per f, è anche possibile lavorare con una base B-spline non vincolata, poiché la penalità in (3) impone automaticamente i vincoli di linearità ai nodi x(1) e x(n) (vedi , Capitolo 5, e , Capitolo 2). Per quanto riguarda la base B-spline, i risultati della stima non dipendono dalla scelta dei nodi di confine: è possibile utilizzare x(1) e x(n) come nodi di confine o includere x(1) e x(n) nell’insieme dei nodi interni.
Se n è grande e l’intervallo è coperto densamente dai dati osservati, di solito non è necessario porre un nodo ad ogni xi,i=1,…,n. Invece, la spline di lisciatura può essere approssimata da una spline di regressione penalizzata che usa un insieme ridotto di nodi. Una classe molto popolare di spline di regressione penalizzate sono le P-spline, che si basano sulla base cubica B-spline e su un insieme ‘grande’ di nodi equidistanti (di solito, 10-40). Invece di valutare l’integrale in (3), le P-spline si basano su una penalità di secondo ordine definita da
che, nel caso di nodi uniformemente spaziati, può essere dimostrato essere un’approssimazione a Jβ. L’operatore di differenza del secondo ordine Δ2 è definito da Δ2βk:=(βk-βk-1)-(βk-1-βk-2). La penalità può quindi essere espressa come β⊤Pβ, dove P è definito da D⊤D con D una matrice di differenze. Si deduce facilmente che lo stimatore risultante di β ha la stessa struttura di 2, con Ω sostituito da P.
Una proprietà conveniente delle P-spline è che sono numericamente stabili e molto facili da definire e implementare. In particolare, è molto più facile impostare la matrice di differenza D che la matrice Ω. Inoltre, è semplice estendere la penalità Jβ (e quindi la matrice D) alle differenze di ordine superiore Δq con q>2. È anche possibile utilizzare una sequenza di nodi che non è uniformemente spaziata; in questo caso, è necessario introdurre dei pesi. Poiché le P-spline con nodi non uniformemente spaziati sono raramente usate nella pratica, non le consideriamo qui e ci riferiamo invece a
Le spline liscianti e le P-spline superano il problema della selezione dei nodi in qualche misura. La loro filosofia è di usare un gran numero di nodi e poi lasciare che λ controlli la quantità di morbidezza. Questo si traduce in un ulteriore parametro di regolazione, senza un consenso generale su come regolare questo parametro. Alcuni modi popolari per determinare il valore “ottimale” di λ usano la convalida incrociata generalizzata (GCV), l’AIC o una rappresentazione a modello misto.
Splines in R
Il pacchetto di installazione base di R contiene un insieme di funzioni che possono adattare semplici spline polinomiali e spline di smussamento. Ulteriori funzioni sono incluse nella libreria splines scritta da DM Bates e WN Venables. Il pacchetto è stato il cavallo di battaglia dell’adattamento di spline per molti anni e ora fa parte della distribuzione base di R. Ci sono più di 100 altri pacchetti che dipendono da splines quando vengono caricati. Il pacchetto contiene diverse funzioni per creare basi di spline, come bs per le B-spline e ns per le spline naturali, che sono ampiamente utilizzate, ma anche alcune funzioni più specializzate per creare funzioni di base (come periodicSpline che crea una spline di interpolazione periodica) o comandi che sono utili come il comando predict.bSpline che valuterebbe una spline a nuovi valori di X.
I valori predefiniti di bs creeranno una base B-spline cubica con due nodi di confine e un nodo interno posto alla mediana dei valori dei dati osservati. Una maggiore flessibilità può essere ottenuta dall’utente, aumentando il posizionamento e il numero di nodi e/o cambiando la loro posizione. La figura 3 (codice 2 nel file aggiuntivo 1: Appendice) mostra le B-spline create con diverse opzioni. La parte superiore presenta spline lineari, cioè polinomi del primo ordine (il grado è uno) collegati tra loro su nodi equidistanti. La parte inferiore presenta polinomi cubici (grado 3).
Si noti che le B-spline create in R con bs() sono automaticamente delimitate dall’intervallo dei dati, e che i nodi aggiuntivi (τ1,…,τd) sono impostati uguali ai nodi di confine, dando più nodi alle due estremità del dominio. Questo approccio è utile nei casi univariati e ha alcune caratteristiche computazionalmente interessanti. Tuttavia, se si lavora su un problema di lisciatura bidimensionale, usando prodotti tensoriali di B-spline, o quando si lavora con P-spline, questa base non è adatta e può portare a risultati spuri.
Le spline naturali possono essere create all’interno del pacchetto splines, usando il comando ns. Per default, a meno che l’utente non specifichi i gradi di libertà o i nodi, la funzione restituisce una linea retta all’interno dei nodi di confine. La figura 4 (codice 3 nel file aggiuntivo 1: Appendice) mostra spline naturali create con diverse opzioni.
Per illustrare come queste funzioni possono essere usate in pratica, consideriamo ancora i dati della Sezione 2.0.1. La figura 5 (creata da (codice 4 nel file aggiuntivo 1: Appendice)) mostra i fits ottenuti utilizzando i seguenti comandi: poly() per semplici spline polinomiche ortogonali, smooth.spline() per spline di smoothing, bs() e ns() da library splines, rispettivamente per B-spline e natural splines. Il grafico in alto a sinistra mostra un semplice adattamento lineare ai dati (linea tratteggiata) e un adattamento polinomiale di terzo grado che è in grado di catturare la relazione più complessa tra le variabili. Il grafico in alto a destra è però particolarmente interessante, poiché presenta i fit utilizzando i valori di default delle funzioni spline. La linea verde proviene dalle funzioni poly() e ns() che di default definiscono entrambe una linea retta. All’altro estremo, la linea blu è un adattamento dalla funzione smooth.spline() che se non sono specificati i gradi di libertà tende a smussare i dati, cioè a produrre un adattamento molto flessibile e sinuoso basato -qui- su 45 gradi di libertà. Un adattamento visivamente ragionevole ai dati può essere ottenuto quando vengono specificati quattro gradi di libertà (grafico in basso a sinistra). Si può vedere che ci sono alcune differenze a seconda della base scelta. La base polinomiale (linea nera) è un po’ più flessibile delle altre, specialmente alle età più alte. D’altra parte, una spline di lisciatura limitata a soli quattro gradi di libertà è più rigida di altri approcci, ma probabilmente lima troppo i dati a piccole età, tra gli anni 0 e 10. Tra i due estremi, le B-spline e le spline naturali forniscono un adattamento molto simile che cattura l’effetto delle piccole età e tende ad essere meno influenzato dai casi estremi alla fine dello spettro di età. Infine, il grafico in basso a destra mostra quanto più flessibile diventi l’adattamento con ulteriori gradi di libertà e suggerisce un potenziale errore di adattamento dovuto all’uso di gradi di libertà eccessivi.
Una nota sui gradi di libertà
In pratica, è sempre utile definire una spline per gradi di libertà. Questo approccio è particolarmente utile quando si lavora con le B-spline e le spline naturali. Le B-spline hanno d+K, mentre una funzione base spline cubica naturale con K nodi ha K+1 gradi di libertà, rispettivamente. Per impostazione predefinita, la funzione bs in R crea B-spline di grado 3 senza nodi interni e nodi di confine definiti nell’intervallo della variabile X. Come tale la funzione crea tre funzioni base. Consideriamo ora il seguente caso: quando un utente definisce una B-spline con un nodo interno alla mediana di X (bs(x,nodi=mediana(x)) il software creerà quattro funzioni (d=3 più K=1 nodi interni, quattro gradi di libertà). Se invece l’utente specifica nella funzione i nodi di confine all’interno dell’argomento nodi (bs(x,knots=c(min(x),median(x),max(x))))), la funzione avrà sei gradi di libertà (d =3 più k =3). Una cautela simile dovrebbe essere presa con la funzione ns.
Quando si lavora con le spline di lisciatura, non è facile specificare i gradi di libertà, poiché essi variano a seconda della dimensione della penalizzazione. Tuttavia, in pratica, le spline penalizzate possono anche essere limitate a un numero massimo di gradi di libertà o di gradi di libertà desiderati.
Altri pacchetti spline
In generale, l’elenco esteso dei pacchetti spline contiene sia approcci abbastanza simili a quello presentato qui sia casi molto specializzati che mirano ad applicazioni specifiche. Nella tabella 1 alcuni di questi pacchetti sono presentati insieme al numero di download. Il numero si riferisce al numero di volte che un pacchetto è stato scaricato ma non agli utenti unici. Va oltre lo scopo di questo lavoro descrivere in dettaglio tutti questi approcci.