venerdì 31 luglio 2009

Come eliminare il fattore stagionalità da un'analisi del trend

"Anyone who tries to analyse a time series without plotting it first is asking for trouble" (Chatfield 1996).

L'analisi di un trend, con misurazioni stagionale, o generalmente cicliche, non può essere misurata con una semplice regressione lineare sui dati di cui disponiamo, se non siamo certi che questi dati siano indipendenti dal periodo di tempo in cui sono stati registrati. Per fittare un modello più preciso occorre quindi tenere presente anche dell'importante dato della stagionalità, e aggiustare il modello di conseguenza (si parla di seasonal adjustment.
Osserviamo ad esempio i dati seguenti, che si riferiscono alle misurazioni effettuate nelle 4 stagioni per 5 anni successivi (totale 5x4 = 20 osservazioni):

21, 16, 11, 19, 25, 19, 14, 22, 29, 23, 20, 26, 33, 27, 22, 30, 36, 29, 23, 33

Ovviamente non possiamo continuare la nostra disamina senza seguire il consiglio di Chatfield (si sa, l'occhio umano vede molto più in un grafico di qualsiasi computer al mondo!):


x <- c(1:20)
y <- c(21, 16, 11, 19, 25, 19, 14, 22, 29, 23, 20, 26, 33, 27, 22, 30, 36, 29, 23, 33)

plot(x, y, type="l")





E' evidente come l'andamento dei dati dipenda dalla stagione. Fittiamo il modello e disegnamo la retta di regressione:


fit <- lm(y ~ x)
abline(fit, col="red", lwd=2)
summary(fit)

Call:
lm(formula = y ~ x)

Residuals:
Min 1Q Median 3Q Max
-7.6361 -2.0744 0.1662 2.4632 7.1188

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 15.5789 2.2066 7.060 1.39e-06 ***
x 0.7925 0.1842 4.302 0.000429 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.75 on 18 degrees of freedom
Multiple R-squared: 0.507, Adjusted R-squared: 0.4796
F-statistic: 18.51 on 1 and 18 DF, p-value: 0.000429




L'equazione della retta di regressione è: y = 15.5789 + 0.7925x. Sia l'intercetta, sia il coefficiente di regressione b sono significativi. Ma il multiple R-squared è molto basso, solo del 50.7%. Ciò ci obbligherebbe a rifiutare il modello, essendo insufficientemente adeguato.

Ora, creiamo un vettore contenente i valori predetti dalla retta di regressione, che possono essere facilmente calcolati andando a sostituire i valori delle x nella equazione della retta di regressione trovata. In R si può usare comodamente la funzione predict():

val_pred <- predict(fit)


Adesso sottraiamo ai valori reali (vettore y), i valori predetti dalla equazione (vettore val_pred), ottenendo così un vettore contenente i cosiddetti valori de-trended:

de_trended <- y - val_pred


Adesso riarrangiamo questi valori al fine di ottenere una tabella a doppia entrata: anno e valore de-trended per ogni quarto. Notate che il 1° quarto, dei 5 anni, è contenuto nelle caselle 1, 5, 9, 13, 17; il 2° quarto nelle caselle 2, 6, 10, 14, 18; e così via.

stag <- matrix(de_trended,ncol=4,byrow=TRUE)
colnames(stag) <- c("Quarto1","Quarto2","Quarto3","Quarto4")
rownames(stag) <- c("Anno1","Anno2","Anno3","Anno4","Anno5")
stag <- as.table(stag)
stag

Quarto1 Quarto2 Quarto3 Quarto4
Anno1 4.6285714 -1.1639098 -6.9563910 0.2511278
Anno2 5.4586466 -1.3338346 -7.1263158 0.0812030
Anno3 6.2887218 -0.5037594 -4.2962406 0.9112782
Anno4 7.1187970 0.3263158 -5.4661654 1.7413534
Anno5 6.9488722 -0.8436090 -7.6360902 1.5714286



Per calcolare la media di ogni colonna (ossia la media delle misurazioni per ogni stagione), procediamo così:

mediacol <- colMeans(stag)
mediacol

Quarto1 Quarto2 Quarto3 Quarto4
6.0887218 -0.7037594 -6.2962406 0.9112782


Adesso trasformiamo in vettore questi dati, e poi creiamo un vettore contenente le medie, ripetuto 5 volte (quanto il numero di anni):


mediacolvec <- as.vector(mediacol)
stag_means <- rep(mediacolvec, 5)


Ora creiamo il vettore differenza tra i dati misurati (vettore y) e le medie stagionali (vettore stag_means):


ValAdj <- y - stag_means


Infine ricaviamo i valori residui, sottraendo ai valori misurati (vettore y) le medie stagionali (vettore stag_means) e i valori predetti dal modello (vettore val_pred):


residuals <- y - stag_means - val_pred


Può essere utile riunire tutti i dati calcolati fino ad ora in un unico dataframe:


data <- data.frame(x, y, val_pred, de_trended, stag_means, ValAdj, residuals)
data

x y val_pred de_trended stag_means ValAdj residuals
1 1 21 16.37143 4.6285714 6.0887218 14.91128 -1.4601504
2 2 16 17.16391 -1.1639098 -0.7037594 16.70376 -0.4601504
3 3 11 17.95639 -6.9563910 -6.2962406 17.29624 -0.6601504
4 4 19 18.74887 0.2511278 0.9112782 18.08872 -0.6601504
5 5 25 19.54135 5.4586466 6.0887218 18.91128 -0.6300752
6 6 19 20.33383 -1.3338346 -0.7037594 19.70376 -0.6300752
7 7 14 21.12632 -7.1263158 -6.2962406 20.29624 -0.8300752
8 8 22 21.91880 0.0812030 0.9112782 21.08872 -0.8300752
9 9 29 22.71128 6.2887218 6.0887218 22.91128 0.2000000
10 10 23 23.50376 -0.5037594 -0.7037594 23.70376 0.2000000
11 11 20 24.29624 -4.2962406 -6.2962406 26.29624 2.0000000
12 12 26 25.08872 0.9112782 0.9112782 25.08872 0.0000000
13 13 33 25.88120 7.1187970 6.0887218 26.91128 1.0300752
14 14 27 26.67368 0.3263158 -0.7037594 27.70376 1.0300752
15 15 22 27.46617 -5.4661654 -6.2962406 28.29624 0.8300752
16 16 30 28.25865 1.7413534 0.9112782 29.08872 0.8300752
17 17 36 29.05113 6.9488722 6.0887218 29.91128 0.8601504
18 18 29 29.84361 -0.8436090 -0.7037594 29.70376 -0.1398496
19 19 23 30.63609 -7.6360902 -6.2962406 29.29624 -1.3398496
20 20 33 31.42857 1.5714286 0.9112782 32.08872 0.6601504



Adesso creiamo il modello con i valori aggiustati tenendo conto della stagionalità (vettore ValAdj):


fitAdj <- lm(ValAdj ~ x)
summary(fitAdj)


Call:
lm(formula = ValAdj ~ x)

Residuals:
Min 1Q Median 3Q Max
-2.01489 -0.34255 -0.07942 0.35628 1.96029

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 14.74507 0.37554 39.26 < 2e-16 ***
x 0.87190 0.03135 27.81 3.05e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.8084 on 18 degrees of freedom
Multiple R-squared: 0.9773, Adjusted R-squared: 0.976
F-statistic: 773.5 on 1 and 18 DF, p-value: 3.048e-16



Osservate il valore di multiple R-squared: questo modello così creato ha un valore molto più elevato, addirittura del 97.73%.
Disegnamo anche la retta di regressione aggiustata, e osserviamo il grafico:


abline(fitAdj, lwd=2, col="blue")
text(max(x)-1.5, max(predict(fit))-2, paste("Without seas adj"), pos=1, col="red")
text(max(x)-3.1, max(predict(fitAdj)), paste("With seas adj"), pos=1, col="blue")





1 Commento:

Anonimo ha detto...

prima di tutto complimenti per il blog, è utilissimo e ben strutturato.
Non ho capito una base del Seasonal adjustment. Se ho una variabile che E' stagionale (per esempio le piogge negli anni, o il ciclo fenologico di una pianta) utilizzare questo metodo è scorretto?

francesco

Posta un commento

Statistiche... del blog!

In questo blog ci sono posts e commenti.

Visualizzazioni totali (dal 01.06.2010)

Follow me on...