Regressione segmentale

Qualche tempo fa un mio cliente mi ha chiesto di analizzare i dati relativi alle visite di un suo sito commerciale. Mi sono stati forniti i numeri corrispondenti alle visite giornaliere ricevute da quel sito per 90 giorni (3 mesi). Queste informazioni sono state prese da Google Analytics, esportate in Excel e quindi passate in R, come spiegato in questo post.
Ecco di seguito i dati:

MESE 1: 51, 49, 44, 31, 32, 34, 36, 42, 31, 40, 27, 22, 27, 33, 45, 47, 43, 26, 21, 36, 46, 46, 31, 29, 23, 36, 48, 188, 148, 116
MESE 2: 32, 164, 136, 180, 184, 148, 124, 148, 136, 136, 192, 156, 132, 168, 220, 120, 116, 164, 220, 144, 120, 132, 124, 116, 120, 180, 104, 184, 160, 132
MESE 3: 100, 192, 132, 136, 192, 140, 156, 203, 196, 350, 273, 203, 322, 217, 126, 315, 294, 238, 189, 245, 161, 273, 385, 315, 287, 294, 350, 259, 168, 294

Ciò che mi era stato chiesto, era un'analisi del trend di visitatori. Il numero di visitatori giornalieri tende ad aumentare? A diminuire?
Per cominciare l'analisi, è sempre utile disegnare uno scatterplot con i valori, per farsi un'idea di cosa dobbiamo andare a cercare, e quali valori ci si può attendere. Eccolo di seguito:

y <- c(51, 49, 44, 31, 32, 34, 36, 42, 31, 40, 27, 22, 27, 33, 45, 47, 43, 26, 21, 36, 46, 46, 31, 29, 23, 36, 48, 188, 148, 116, 32, 164, 136, 180, 184, 148, 124, 148, 136, 136, 192, 156, 132, 168, 220, 120, 116, 164, 220, 144, 120, 132, 124, 116, 120, 180, 104, 184, 160, 132, 100, 192, 132, 136, 192, 140, 156, 203, 196, 350, 273, 203, 322, 217, 126, 315, 294, 238, 189, 245, 161, 273, 385, 315, 287, 294, 350, 259, 168, 294)
x <- c(1:length(y))

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



Possiamo osservare chiaramente che l'andamento non è affatto lineare. Oltre alle variazioni dovute al caso, si possono visualizzare tre periodi con frequenze notevolmente differenti. Effettuare quindi una regressione lineare per l'analisi del trend temporale è poco efficace. E' abbastanza prevedibile che il coefficiente angolare risulterà positivo, e che la retta di regressione sarà abbastanza inclinata verso l'alto. Ma come spiegarsi quegli apparenti 3 periodi?

Questo è un classico caso in cui ci si affida all'analisi della regressione segmentale (segmented regression, o breakpoints regressione, o piecewise regression). La regressione segmentale è un metodo utilizzato per studiare l'andamento di certi dati in periodi temporali che possono differire per un qualsiasi motivo. In questo caso però, inizialmente non ero a conoscenza del perchè ci fossero così forti discrepanze tra i 3 periodi, avendo a disposizione solo i tre dati. Tuttavia la presenza di un grafico così accentuato, mi ha portato comunque a procedere con l'analisi della regressione con punti di rottura.
Sostanzialmente i metodi utilizzati per questa analisi sono gli stessi utilizzati per l'analisi della regressione lineare, di cui abbiamo già parlato. In pratica si divide il campione in segmenti, individuando i punti di rottura (breakpoints), e si fa una regressione lineare sui vari segmenti. Ma dove si trova questo punto di rottura?
Possiamo decidere noi dove posizionarlo: se fossi stato a conoscenza, ad esempio, che il giorno X il mio cliente ha apportato una sostanziale modifica al suo sito, quello sarebbe stato un punto di rottura. Non avendo tali informazioni, è necessario stimare statisticamente in corrispondenza di quali punti l'andamento subisce una forte modifica.

In R esiste una funzione apposita, la funzione breakpoints(), disponibile nella libreria strucchange. Dopo aver scaricato la libreria (clickate su Pacchetti > Installa Pacchetti > Selezione del mirror (Italia) > strucchange), bisogna aprire il pacchetto (così le funzioni incluse potranno essere utilizzate), e applichiamo la funzione breakpoints, la quale richiede in input il modello che vogliamo analizzare. In R quindi si procede così:

library(strucchange)

bp <- breakpoints(y ~ x)
> bp

Optimal 3-segment partition:

Call:
breakpoints.formula(formula = y ~ x)

Breakpoints at observation number:
27 69

Corresponding to breakdates:
0.3 0.7666667

L'output ci dice che possiamo dividere i dati in esame in 3 segmenti (come avevamo previsto). I punti di rottura corrispondono ai valori nelle posizioni 27 e 69. Procederemo quindi con la stima della regressione nei tre intervalli (1,27), (28,69), (70,90).

Innanzitutto disegnamo due linee verticali sul plot precedentemente creato, in corrispondenza dei due breakpoints. La funzione abline fa al caso nostro. Essa aggiunge una linea al grafico già in finestra; specificando il parametro v, si definisce una linea verticale con coordinata x corrispondente; il comando lty=4 disegna una linea tratto-punto (dash-dot):

abline(v=27, lwd=2, lty=4)
abline(v=69, lwd=2, lty=4)

A questo punto creiamo i modelli di regressioni lineare totale, e sui 3 segmenti:

fit_totale <- lm(y ~ x)
fit1 <- lm(y[1:27]~x[1:27])
fit2 <- lm(y[28:69]~x[28:69])
fit3 <- lm(y[70:90]~x[70:90])

Fatto ciò, possiamo analizzare in dettaglio, con la funzione summary(), i 3 modelli, ma in questo caso è altrettanto utile porre in grafico le rette di regressione. Per disegnare la retta di regressione lineare totale, è sufficiente utilizzare ancora una volta la funzione abline:

abline(fit_totale, lwd=2, col="red")

Adesso dobbiamo disegnare solo 3 segmenti, e non le rette di regressione, dei 3 modelli segmentali. Non possiamo pertanto utilizzare la funzione abline, ma useremo la funzione lines. Questa funzione disegna una linea unendo i punti di coordinate (x0, y0) e (x1,y1). Dobbiamo fornire quindi alla funzione due vettori, contenenti i valori (x0, x1) e (y0, y1) (quindi non le coordinate dei due punti, ma prima le X poi le Y), affinchè i due punti vengano uniti. Vediamo come ricavare le coordinate.

Abbiamo detto che il modello fit1 è valido nel tratto delle x (che sono i giorni) tra 1 e 27. L'equazione generale di una retta è y = a + bx. La variabile indipendente x assume i valori 1 (x0) e 27 (x1); calcoliamo i valori y0 e y1.
In un qualsiasi modello fittato, abbiamo che:
a = coef(modello)[1]
b = coef(modello)[2]

Quindi in generale y0 sarà uguale a a + b * x0, e y1 sarà pari a a + b * x1. Tradotto in R tutto questo diventa:

x01 <- 1
x11 <- 27
y01 <- coef(fit1)[1] + coef(fit1)[2]*x01
y11 <- coef(fit1)[1] + coef(fit1)[2]*x11
x_var1 <- c(x01, x11)
y_var1 <- c(y01, y11)

Per il secondo modello, abbiamo:

x02 <- 27
x12 <- 69
y02 <- coef(fit2)[1] + coef(fit2)[2]*x02
y12 <- coef(fit2)[1] + coef(fit2)[2]*x12
x_var2 <- c(x02, x12)
y_var2 <- c(y02, y12)

E infine per il terzo modello:

x03 <- 69
x13 <- 90
y03 <- coef(fit3)[1] + coef(fit3)[2]*x03
y13 <- coef(fit3)[1] + coef(fit3)[2]*x13
x_var3 <- c(x03, x13)
y_var3 <- c(y03, y13)

Adesso siamo pronti ad utilizzare la funzione lines() inserendo le coordinate:

lines(x_var1, y_var1, lwd=3, col="blue")
lines(x_var2, y_var2, lwd=3, col="blue")
lines(x_var3, y_var3, lwd=3, col="blue")

Ecco come appare il grafico completo:



Osserviamo la linea rossa (retta di regressione totale): mostra un andamento crescente notevolmente marcato. Ma ora osserviamo i tre segmenti di regressione: il primo mostra un andamento discendente, il secondo appena crescente, il terzo un pò più crescente. Quindi la retta di regressione totale avrebbe fuorviato la nostra analisi: sembrerebbe che il sito sia in crescita continua, mentre così non è.
Come spiegare questa diversità in termini di risultati? Interrogando il mio cliente, ho scoperto che proprio nei giorni che la funzione breakpoints aveva predetto (il giorno 27 e il giorno 69), egli aveva apportato grosse modifiche al proprio sito, e ne aveva informato i propri clienti fedeli tramite newsletter. Questo ha fatto notevolmente aumentare il numero di visitatori (grazie al passaparola), ma questi rimangono pressochè costanti, fino al giorno in cui una nuova modifica non viene effettuata. Allora si ha un nuovo giro di parole, nuovi visitatori, e le visite aumentano di nuovo.

Il mio consiglio finale al cliente?
Apportare ogni giorno modifiche al proprio sito :-D (ovviamente scherzavo!!)

Per completezza logicamente bisogna valutare la significatività dei tre modelli di regressione (con la funzione summary(), e poi da lì si analizza la significatività dei coefficienti). E' inoltre opportuno calcolare gli intervalli di confidenza per i breakpoints stimati, con la funzione confint(bp, breaks=2).
Per una regressione segmentale può essere utile la funzione segmented() presente nella library segmented, che va scaricata come visto sopra. Questa funzione effettua una regressione segmentale, così come abbiamo visto in questo post, stimando i breakpoints.