LOESS

Il metodo dei polinomi locali noto anche come metodo LOESS, abbreviazione di LOcal regrESSion, viene impiegato per adattare una superficie di regressione ai dati attraverso una procedura di smoothing (lisciamento).

L’analisi di regressione è un modo per quantificare la relazione esistente tra una variabile di risposta Y ed una o più variabili esplicative. Si supponga che, data una variabile di risposta Y ed un insieme di \(p\) variabili esplicative \(x = (x_1, …, x_p)\) fissate, la relazione fra \(Y\) e \(x\) sia data da: \[Y = m(x) + \epsilon\] Nel contesto dei trend spaziali x è da intendersi come il vettore s delle coordinate di un punto nella regione di studio D, mentre \(\epsilon\) costituisce il termine di errore. L’assunzione di una struttura polinomiale del predittore può essere eccessivamente semplificatrice se mantenuta fissa su tutto il modello mentre in LOESS questa viene considerata localmente su un intorno di x. In un primo momento si spiegherà questo metodo nel caso univariato per poi estendere la definizione al caso multivariato (che difatti è quello spaziale). Il metodo LOESS permette di ottenere un valore \(\hat y_0\) differenziando il ruolo delle rilevazioni presenti nel campione e in modo tale che le osservazioni più vicine a \(x_0\) abbiano un peso maggiore da quelle che invece si trovano da esso distanti. In questo approccio occorre quindi definire un’opportuna metrica delle esplicative per quantificare la distanza tra \(x_0\) el le altre \(x_k\) con \(k= 1,…,n\). Si supponga che la relazione tra il valore atteso della risposta e c sia data in un intorno ragionevolmente piccolo di z da un polinomio di grado \(q\), quindi il modello viene formalizzato come:

\[m(x) = \beta_0 + \beta_1(x) \space x + \beta_2(x) \space x^2 + … + \beta_q(x) \space x^q\]

Per calcolare la stima di $m (x) $ occorre calcolare le stime dei parametri del polinomio \(\beta_j(x)\). Qeste stime sono ottenute applicando i minimi quadrati ponderati e, quindi, i \(\hat \beta_j(x)\) sono i valori \(\beta_j =\beta_j(x)\) che minimizzano: \[ \sum \bigg\{w_k(x) \space \space (y_k -\beta_0 -\beta_1 x_k - … - \beta_q x_k^q)^2 \bigg\}\] Sove I valori $ w_k(x)$ sono degli opportuni pesi definiti in basse alla distanza di ciascun punto del campione \(x\)

Per questo si ricorre ad opportune funzioni nucleo dette anche, utilizzando la denominazione inglese ormai generalmente diffusa, kernel. La funzione kernel deve avere alcuni requisiti che la rendano idonea a definire un sistema di ponderazione. In particolare, deve essere positiva, con un massimo in 0 e simmetrica rispetto a questo valore. In via di principio qualsiasi funzione con queste caratteristiche potrebbe essere utilizzata per definire i pesi. Tuttavia, la funzione dei pesi cui tipicamente si fa riferimento nel LOESS è la funzione tricubica:

\[W(u) = \begin{cases} (1-|u|^3)^3, & \mbox{se } -1 \leq u \leq 1 \\ 0, & \mbox{altrimenti} \end{cases}\]

Questa funzione è una funzione a supporto compatto, con un massimo in 0 e simmetrica rispetto a questo valore.

Il peso da assegnare alla k-esima osservazione campionaria nella stima della regressione in x è ottenuto come : \[w_k(x) = W\bigg(\frac{x_k-x}{h}\bigg)\]

Dove h è un parametro positivo che rappresenta la finestra di lisciamento, ovvero l’ampiezza dell’intervallo scelto attorno a x entro cui i punti campionari assumono un peso positivo nella stima di m(x). Nella stima LOESS, diversamente da altri metodi di lisciamento, la scelta della finestra non è fissa, ma varia a seconda del punto x. Più in particolare, h viene scelto fissando la percentuale di punti \(\frac{c}{n}\) (detta span) utilizzati per il lisciamento, in modo tale che l’ampiezza della finestra si modifica da x a x al fine di garantire che venga utilizzato sempre lo stesso numero c di valori campionari nella stima di m(x) per diverse scelte di x. Tanto più lo span è piccolo, tanto più locale sarà la stima e, quindi, meno liscia la superficie. Tanto più grande è il valore dello span tanto più liscia sarà la superficie prodotta e tanto “meno locale” la stima.

Se x è un vettore in \(R^p\) e \(d(u,v)\) una metrica di ale spazio per esempio la distanza euclidea tra due punti allora il peso di \(w_k(x)\) da attribuirsi all’osservazione (y_k,x_k) ciene definita tramite la funzione \(w_k(x)=W\bigg(\frac{d(x,x_k)}{h}\bigg)\)

Altre funzioni che si possono utilizzare come kernel

Variogramma empirico (Introduzione)

É possibile capire meglio il Variogramma guardando la nuvola di punti che si genera rappresentando graficamente le semivarianze di tutte le coppie di osservazioni:

Ogni punto rapprensenta un paio di dati osservati, sull’asse X c’è la distanza (metri), mentre sull’asse Y c’è la varianza. La nube di punti contiene tutte le relazioni spaziali nei dati per tutte le possibili distanze tra i campioni, ma non è una funzione continua. E’ difficie interpretarla e comprendere l’esistenza di correlazioni spaziali. Si rende necessario prendere in considerazione un numero ristretto di lag (distanze). In questo modo sono più riconoscibili gli outlier ed è possibile modellizzare la distribuzione spaziale

Il variogramma è semplicemente una rappresentazionesintetica (media) della nuvola di punti delle semivarianze vs distanze. Una funzione per modellare è necessaria per ottenereun valore di varianza per ogni possibile distanza tra i dati alle posizioni note e le posizioni che devono essere predette.

• Il nugget rappresenta le variazioni casuali, della variabilità su piccola scala (inferiore cioè alla scala del campionamento), la variazione che ci sarebbe tra due punti più vicini rispetto alla distanza minima dei punti campionati.
• Il range corrisponde alla distanza in cui il semivariogramma (in un modello stazionario) raggiunge un valore costante detto sill. Intuitivamente è la massima distanza per cui due punti possono essere considerati spazialmente correlati.

• La soglia (sill) rappresenta la varianza massima tra i punti misurati, che per modelli stazionari corrisponde alla varianza di campionamento.

  
  
  

  Come si evidenzierà meglio in seguito il variogramma empirico è uno stimatore del variogramma del processo che risulta adeguato in situazioni di stazionarietà intrinseca

  Effetto pepita (Nugget effect)

  \(\hat \gamma(h)\) assume valori diversi da 0 in prossimità di h = 0 pur essendo per definizione $(0)=0 $ quindi anche a piccole distanze ci sono forti variazioni nella variabile rilevata (isotropia).

  

  Hole effect

  caduta rapida dei valori del variogramma empirico quando la distanza aumenta valori a maggiori distanze sono loro più simili di valori a distanze ridotte, possibili cause sono: presenza di trend, presenza di valori anomali, correlazioni negative o periodicità.

  

  

  La stima del variogramma può essere realizzata seguendo un approccio di tipo non parametrico o parametrico. Il parte seguente è dedicata alla stima non parametrica del variogramma e la successiva alla stima parametrica.

Il metodo dei momenti (MOM) [Matheron, 1962]

I variogramma empirico mette in relazione i valori osservati delle variabili di interesse con la posizione reciproca delle locazioni in cui sono rilevate. In particolare confronta per differenza coppie di valori presi a distanza h:

\[2 \hat \gamma (h) = \frac{1}{\#N(h)} \space \space \sum_{s_i,s_j \in N(h)}\bigg(Y(s_i)-Y(s_j)\bigg)^2 \]

… e di conseguenza il semi-variogramma risulta:

\[ \hat \gamma (h) = \frac{1}{2\space \space \#N(h)} \space \space \sum_{s_i,s_j \in N(h)}\bigg(Y(s_i)-Y(s_j)\bigg)^2 \] dove \(N(h)\) contiene tutte le locazioni che distano h tra loro

Proprietà:

• \(\hat \gamma(h)=\hat \gamma (-h)\) (come per il variogramma teorico)
• Correttezza:\[ E[\hat \gamma (h)]= \hat \gamma (h)\]
• Consistenza:\[Corr\bigg\{ \bigg(Y(s + h)-Y(s)\bigg)^2 \space \space , \space \space \bigg(Y(s' + h)-Y(s')\bigg)^2\bigg\} \approx 0\] Da questa proprietà di consistenza ne deriva che per la legge dei grandi numeri \[\hat \gamma (h) \rightarrow^d \space \space Normale \space Multivariata\] Questa convergenza è tanto piu veloce quanto la proprietà di consistenza è verificata, ossia, avviene lentamente più le osservazioni sono correlate mentre se le osservazioni risultano incorrelate allora il semivariogramma convergerà in distribuzione ad una Normale Multivariata in modo più rapido

In R si utilizza la funzione variog nel pacchetto geoR inserendo il dataset trasformato in un elemento geodata e imponento il parametro estimator.type = “classical”

Stima robusta del varogramma empirico [Cressie e Hawkins, 1980]

I variogramma empirico mette in relazione i valori osservati delle variabili di interesse con la posizione reciproca delle locazioni in cui sono rilevate. In particolare confronta per differenza coppie di valori presi a distanza h:

\[2 \gamma^* (h) = \bigg[ \frac{1}{\#N(h)} \space \space \sum_{s_i,s_j \in N(h)}\bigg|Y(s_i)-Y(s_j)\bigg|^ {\frac 1 2}\bigg]^4 \]

Questo metodo nasce con l’obiettivo di ridurre l’impatto degli outlier. La potenza \(4^a\) consente di ritornare alla scala di \(2 \gamma^* (h)\).

Per ottenere una versione più resistente si utilizza al posto della media la mediana:

\[2 \tilde \gamma (h) = \frac{Med \bigg[ \bigg|Y(s_i)-Y(s_j)\bigg|^{\frac 1 2}\bigg]^4}{B(h)} \\ dove: \\ B(h)= 0.457 \rightarrow \mbox{fattore di correzione} \\ s_i,s_j \in N(h)\]

In R si utilizza la funzione variog nel pacchetto geoR inserendo il dataset trasformato in un elemento geodata e imponento il parametro estimator.type = “modulus”, in questa funzione lo stimatore risulta leggermente lo stesso ma viene calcolato utilizzando le proprietà del Modified Z-Scores:

$$\mathrm{\gamma <!--\; \gamma \; -->}(h)=\frac{[\frac{1}{N_h}\underset{i=1}{\overset{N_h}{\sum <!--\; \sum \; -->}}|Y(x_{i+h})-<!--\; -\; -->Y(x_i){|}^{\frac{1}{2}}{]}^4}{0.914+\frac{0.988}{N_h}}$$

Stimatore kernel del variogramma [Nadaraya-Watson, 1964]

\[2 \bar \gamma (h,u) = \sum_{i, \space j}^{\#N(h)}\space w_{i\space j}\bigg(Y(s_i)-Y(s_j)\bigg)^2 \\ con: \\ \sum w_{i, \space j}=1\\0<w_{i,\space j}<1 \]

I pesi sono ottenuti dalla funzione kernel \(K(u)\)

\[ w_{i, \space j}=w(s_i,s_j,d)= \frac{K\bigg((h- ||s_i - s_j||)/d\bigg)}{\sum_{s_i,s_j \in N(h)} K\bigg((h- ||s_i - s_j||)/d\bigg)} \]

\(d\): è il parametro di lisciamento.

In R si utilizza la funzione variog nel pacchetto geoR inserendo il dataset trasformato in un elemento geodata e imponento il parametro op=“sm”

Stima di massima verosimiglianza del variogramma (ML)

Questo metodo si basa sulla massimizzazione della funzione di verosimiglianza del campione. Occorre quindi specificare la distribuzione congiunta delle osservazioni campionarie. Tipicamente questo approccio è relativamente agevole assumendo che il processo \(Y(s)\) sia un PSS gaussiano (isotropico). Solitamente la distribuzione congiunta del campione è una normale multivariata il cui vettore delle medie e la cui matrice di varianza e covarianza sono data dalla funzione della media e dal covariogramma del processo.
La funzione di log-verosimiglianza tipicamente assume una forma complessa e l’ottimizzazione non può essere risolta per via analitica, ma viene ottenuta numericamente tramite procedure iterative. L’approccio usualmente adottato in questo caso è quella della profilatura della verosimiglianza.

In R si utilizza la funzione likfit nel pacchetto geoR inserendo il dataset trasformato in un elemento geodata e imponento il parametro cov.model=“exponential” e inizializzando i due paramtri (media e covarianza) con ini.cov.pars=c(par1,par2)

IMPLEMENTAZIONE IN R
likfit(dat.om.geo,ini.cov.pars=c(par1,par2),
cov.model=“exponential”,trend=“cte”, fix.nugget=FALSE,nugget=20,
nospatial=FALSE,messages=T)

Stima dei minimi quadrati del variogramma

Il metodo dei minimi quadrati, utilizzato in questo contesto consiste nell’interpolazione del variogramma empirico tramite un modello parametrico. In questo caso non occorre dare alcuna assunzione circa la distribuzione del processo a parte richiedere che sia intrinsecamente stazionario.

Per ottenere lo stimatore in questione serve per prima cosa stimare preliminarmente il variogramma empirico poi serve una valutazione grafica per individuare una o più famiglie parametriche che ne colgano le caratteristiche salienti. infine, si procede all’interpolazione della curva cercando la stima del vettore \(\theta\) che rende minima la distanza tra i punti del variogramma empirico e la curva teorica scelta. piu precisamente calcolato il variogramma empirico a distanze prescelte si considera il vettore dei valori dipendenti dal vettore dei parametri ignoti assunti dal modello parametrico scelto in corrispondenza delle distanze prefissate.

IMPLEMENTAZIONE IN R
 variofit(geodata,
        ini.cov.pars=c(par1,par2),
          cov.model="exponential",
          fix.nugget=FALSE,nugget=20,
          weights = "equal" )

IMPLEMENTAZIONE IN R
 variofit(geoadata,
        ini.cov.pars=c(par1,par2),
          cov.model="exponential",
          fix.nugget=FALSE,nugget=20)