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")




giovedì 30 luglio 2009

Test del chi-quadro per il trend (test di Armitage)

Test di verifica di un trend lineare in proporzioni e frequenze.

Supponiamo di disporre di una tabella di contingenza 2xK, ad esempio relativa al numero di cavie morte dopo somministrazione di dosi crescenti di un certo farmaco sperimentale, oppure relativa alll'incidenza di patologie in gruppi di individui appartenenti a classi di età sempre maggiori. Con il test di Armitage (Cochran-Armitage test for trend; in alcuni software come SPSS è chiamato anche Mantel-Haenszel linear-by-linear association chi-squared test) possiamo verificare se è presente o meno un trend lineare delle proporzioni. L'ipotesi nulla è la presenza di un trend lineare, contro l'ipotesi alternativa di allontanamento da esso. Vediamo come procedere analizzando la seguente tabella:



Per prima cosa dobbiamo verificare se le proporzioni sono tra loro simili, con il test del Chi-quadro per l'indipendenza delle variabili, già presentato in questo post.

morti <- c(19, 29, 24, 20)
tot <- c(516, 589, 293, 167)
prop.test(morti, tot)

4-sample test for equality of proportions without continuity correction

data: con out of tot
X-squared = 19.5233, df = 3, p-value = 0.0002131
alternative hypothesis: two.sided
sample estimates:
    prop 1     prop 2     prop 3     prop 4
0.03682171 0.04923599 0.08191126 0.11976048


Il Chi-quadro ottenuto è pari a 19.5233, e p-value < 0.05. Quindi rifiutiamo l'ipotesi nulla di uguaglianza delle proporzioni. Le 4 proporzioni sono statisticamente differenti.

Procediamo con un test di correlazione per variabili qualitative, con il test del chi-quadro di Pearson illustrato in questo post. Procediamo in questo modo:

tabella = matrix(c(19, 29, 24, 20, 497, 580, 269, 147), nrow=2, byrow=TRUE)
chisq.test(tabella)

        Pearson's Chi-squared test

data:  tabella
X-squared = 20.1478, df = 3, p-value = 0.0001582


Il risultato di questo test ci dice che esiste una correlazione tra le variabili. Quindi è lecito supporre un andamento lineare di queste, e perciò esaminiamo la presenza di un trend lineare con il test di Armitage:

morti <- c(19, 29, 24, 20)
tot <- c(516, 589, 293, 167)
prop.trend.test(morti, tot)

        Chi-squared Test for Trend in Proportions

data:  morti out of tot ,
using scores: 1 2 3 4
X-squared = 18.2109, df = 1, p-value = 1.977e-05


Abbiamo ottenuto un chi-quadro elevato, e p-value < 0.05. Rifiuto l'ipotesi nulla, e quindi posso affermare che esiste un trend lineare significativo delle proporzioni.
Abbiamo quindi individuato che la componente lineare spiega l'andamento delle proporzioni. Ma è questa l'unica spiegazione all'andamento delle proporzioni, o esiste un'altra regola sottesa? Analizziamo la deviazione dalla linearità, calcolando il chi-quadro residuo, che è pari al chi-quadro totale (cioè il chi-quadro di Pearson), meno il chi-quadro lineare (cioè il chi-quadro di Armitage):

corr <- chisq.test(tabella)
trend <- prop.trend.test(morti, tot)

chi_quadro_residuo <- corr$statistic - trend$statistic

> chi_quadro_residuo

X-squared
1.936932


Il chi-quadro residuo ha gradi di libertà pari a k-2, quindi 2 in questo caso. Calcoliamo quindi il chi-quadro tabulato:

> qchisq(0.950, 2)
[1] 5.991465


Poichè il chi-quadro residuo è inferiore al valore tabulato, accetto l'ipotesi nulla: la componente residua non è significativa, quindi la componente lineare è la sola a spiegare l'andamento delle proporzioni (se avessi rifiutato, significherebbe che la componente lineare non è l'unica variabile esplicativa, ma ne esistono anche altre da ricercare).

Grafici con R: diagramma circolare (o grafico a torta)

Questo è forse uno dei grafici più comuni, essendo una rappresentazione grafica di veloce lettura anche da parte dei meno attenti.
Vediamo alcuni modi per rappresentare un grafico a torta in R. Utilizzerò qui i dati relativi alle elezioni dell'EuroParlamento del 2009: su un totale di 736 europarlamentari, questi i numeri dei singoli partiti presenti:

EPP: 265
S&D: 184
ALDE: 84
Greens-EFA: 55
ECR: 54
EUL-NGL: 35
EFD: 32
NI: 27


Creiamo un grafico a torta semplice con i seguenti comandi:

slices <- c(265, 184, 84, 55, 54, 35, 32, 27)
lbls <- c("EPP", "S&D", "ALDE", "Greens-EFA", "ECR", "EUL-NGL", "NFD", "NI")
pie(slices, labels = lbls, main="Grafico a torta dei partiti presenti all'EuroParlamento 2009")




Come vedete dall'immagine, sono già presenti di default le etichette per ogni fetta del grafico. Può essere utile aggiungere le percentuali:

slices <- c(265, 184, 84, 55, 54, 35, 32, 27)
lbls <- c("EPP", "S&D", "ALDE", "Greens-EFA", "ECR", "EUL-NGL", "NFD", "NI")
pct <- round(slices/sum(slices)*100) #calcolo delle percentuali
lbls <- paste(lbls, pct) # aggiungo il numero percentuale alle etichette
lbls <- paste(lbls,"%",sep="") # aggiungo il simbolo % alle etichette
pie(slices, labels = lbls, main="Grafico a torta dei partiti presenti all'EuroParlamento 2009")




Aggiungiamo, tra parentesi, anche il numero netto di EuroParlamentari per ogni partito:

slices <- c(265, 184, 84, 55, 54, 35, 32, 27)
lbls <- c("EPP", "S&D", "ALDE", "Greens-EFA", "ECR", "EUL-NGL", "NFD", "NI")
pct <- round(slices/sum(slices)*100) #calcolo delle percentuali
lbls <- paste(lbls, pct) # aggiungo il numero percentuale alle etichette
lbls <- paste(lbls,"%",sep="") # aggiungo il simbolo % alle etichette
lbls <- paste(lbls, "(", slices, ")") #aggiungo tra parentesi il numero di europarlamentari
pie(slices, labels = lbls, main="Grafico a torta dei partiti presenti all'EuroParlamento 2009")




Senza ulteriori specificazioni, ogni fetta della torta viene colorata con un colore leggermente diverso. Possiamo cambiare questa impostazione, ad esempio scegliendo un unico colore per ogni fetta (è sufficiente aggiungere, ad esempio, col="red"), oppure possiamo decidere di colorare diversamente ogni fetta, con il seguente comando: col=rainbow(length(lbls)). Ecco come appare ora il grafico:



Oppure possiamo decidere noi singolarmente il colore di ogni fetta. Ad esempio se vogliamo far alternare i colori rosso e blu, creiamo un vettore contenente i colori che vogliamo utilizzare, nell'ordine delle fette:

color <- c("red", "blue", "red", "blue", "red", "blue", "red", "blue")
pie(slices, labels = lbls, col=color, main="Grafico a torta dei partiti presenti all'EuroParlamento 2009")


Se vogliamo invertire l'ordine con cui sono presentati i dati nella torta (normalmente in senso antiorario), aggiungiamo clockwise=TRUE tra gli attributi di pie.




Per creare un grafico a torta tridimensionale, dobbiamo caricare il package plotrix, ed utilizzare la funzione pie3D(), con gli stessi attributi visti finora, ed altri:

library(plotrix)
slices <- c(265, 184, 84, 55, 54, 35, 32, 27)
lbls <- c("EPP", "S&D", "ALDE", "Greens-EFA", "ECR", "EUL-NGL", "NFD", "NI")
pct <- round(slices/sum(slices)*100) #calcolo delle percentuali
lbls <- paste(lbls, pct) # aggiungo il numero percentuale alle etichette
lbls <- paste(lbls,"%",sep="") # aggiungo il simbolo % alle etichette
pie3D(slices, labels = lbls, main="Grafico a torta dei partiti presenti all'EuroParlamento 2009")




Questo grafico può essere ampiamente personalizzato. Possiamo ad esempio specificare l'altezza della torta (height=) o il suo raggio (radius=), l'angolo a partire dal quale cominciare a tagliare le fette espresso in radianti (start=), a quale angolatura presentare la torta (theta=), i colori (col=), la posizione delle etichette (labelpos=), e il loro colore (labelcol=), e possiamo scegliere di separare le fette con il comando explode=. Esaminate questo codice di esempio in R:

pie3D(slices, labels = lbls, height=.1, radius=.9, explode=.2, labelcol="black", main="Grafico a torta dei partiti presenti all'EuroParlamento 2009", theta=pi/12)

martedì 28 luglio 2009

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.

domenica 26 luglio 2009

Tecniche di regressione polinomiale

Nel software statistico R esistono numerosi pacchetti e numerose funzioni per stimare un modello di regressione polinomiale lineare. Questo è però uno degli argomenti più oscuri, a mio avviso, nelle guide e nei tutorial di R, i quali trattano solo di sfuggita questo capitolo. Un altro problema che ho dovuto affrontare, con non poche difficoltà, è la rappresentazione grafica dei modelli trovati. Mi duole ammettere che R (il software statistico che personalmente preferisco) ha non poche limitazioni per quanto riguarda questo specifico argomento.
Vi chiedo scusa per questo sfogo personale, ma ho perso quasi un pomeriggio per giungere alle conclusioni che ora vi presenterò :-)

Supponiamo di voler creare un polinomio che possa approssimare al meglio il seguente dataset, relativo alla popolazione di una certa città italiana nell'arco di 10 anni. In tabella sono riassunti i dati:



Per prima cosa importiamo i dati in R. Se questa è una tabella di Excel, possiamo utilizzare la funzione read.delim("clipboard"), come discusso in un post precedente. Otterremo così un dataframe contenente tutti i valori (questa è la scelta che io consiglio, cioè quella di creare un dataframe). In alternativa possiamo inserire i dati nel modo consueto:

Anno <- c(1959, 1960, 1961, 1962, 1963, 1964, 1965, 1966, 1967, 1968, 1969)
Popolazione <- c(4835, 4970, 5085, 5160, 5310, 5260, 5235, 5255, 5235, 5210, 5175)


A questo punto creiamo il dataframe di nome data:

data <- data.frame(Anno, Popolazione)
data
   Anno Popolazione
1  1959        4835
2  1960        4970
3  1961        5085
4  1962        5160
5  1963        5310
6  1964        5260
7  1965        5235
8  1966        5255
9  1967        5235
10 1968        5210
11 1969        5175


Io consiglio di procedere lavorando con i dataframe, ma questa scelta non è obbligatoria: potremmo anche lavorare con i vettori.

A questo punto può essere utile mettere in grafico questi valori, per osservare l'andamento e farci un'idea della forma del grafico finale. Per comodità modifichiamo la colonna Anno che contiene numeri molto elevati, creando un intorno dello zero, in questo modo:

data$Anno <- data$Anno - 1964
data
   Anno Popolazione
1    -5        4835
2    -4        4970
3    -3        5085
4    -2        5160
5    -1        5310
6     0        5260
7     1        5235
8     2        5255
9     3        5235
10    4        5210
11    5        5175


Richiamiamo il grafico così:

plot(data$Anno, data$Popolazione, type="b")



A questo punto possiamo cominciare con la ricerca di un modello polinomiale che approssimi adeguatamente i nostri dati. Innanzitutto specifichiamo che quello che vogliamo è un polinomio in funzione di X, ossia un raw polynomial in inglese, che è differente dall'orthogonal polynomial. Questa è una precisazione importante perchè i comandi e i risultati cambieranno nei due casi in R. Quindi noi cerchiamo una funzione di X del tipo:



A che grado del polinomio fermarci? Dipende dal grado di precisione che cerchiamo. Maggiore è il grado del polinomio, maggiore sarà la precisione del modello, ma maggiori saranno le difficoltà di calcolo; inoltre bisogna verificare la significatività dei coefficienti che vengono trovati. Ma veniamo subito al dunque.

In R per fittare un modello di regressione polinomiale (non ortogonale) esistono due metodi, tra loro identici. Supponiamo di cercare i valori dei coefficienti beta per un polinomio di 1° grado, quindi di 2° grado, quindi di 3° grado:

fit1 <- lm(data$Popolazione ~ data$Anno)
fit2 <- lm(data$Popolazione ~ data$Anno + I(data$Anno^2))
fit3 <- lm(data$Popolazione ~ data$Anno + I(data$Anno^2) + I(data$Anno^3))


Oppure possiamo scrivere più rapidamente per i polinomi di grado 2 e 3:

fit2b <- lm(data$Popolazione ~ poly(data$Anno, 2, raw=TRUE))
fit3b <- lm(data$Popolazione ~ poly(data$Anno, 3, raw=TRUE))


La funzione poly che vediamo applicata è utile qualora si desideri ottenere un polinomio di grado elevato, perchè evita di scrivere esplicitamente la formula. Se specifichiamo raw=TRUE, i due metodi visti forniscono lo stesso output; se invece non specifichiamo raw=TRUE, la funzione poly ci fornirà i valori dei parametri beta di un polinomio ortogonale, che è differente dalla formula generale che ho scritto sopra, sebbene i modelli siano entrambi efficaci. Rimando a testi di statistica per la differenza tra le due metodiche.

Analizziamo ora l'output.

> summary(fit2)

Call:
lm(formula = data$Popolazione ~ data$Anno + I(data$Anno^2))

Residuals:
    Min      1Q  Median      3Q     Max
-46.888 -18.834  -3.159   2.040  86.748

Coefficients:
               Estimate Std. Error t value Pr(>|t|)
(Intercept)    5263.159     17.655 298.110  < 2e-16 ***
data$Anno        29.318      3.696   7.933 4.64e-05 ***
I(data$Anno^2)  -10.589      1.323  -8.002 4.36e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 38.76 on 8 degrees of freedom
Multiple R-squared: 0.9407, Adjusted R-squared: 0.9259
F-statistic: 63.48 on 2 and 8 DF, p-value: 1.235e-05


L'output di summary(fit2b) è assolutamente identico. Abbiamo ottenuto i valori di beta0 (5263.159), beta1 (29.318) e beta2 (-10.589), i quali risultano essere tutt'e 3 significativi. L'equazione del polinomio di 2° grado del nostro modello è quindi:



Se cercassimo un polinomio di 3° grado, avremmo:

> summary(fit3)

Call:
lm(formula = data$Popolazione ~ data$Anno + I(data$Anno^2) +
    I(data$Anno^3))

Residuals:
    Min      1Q  Median      3Q     Max
-32.774 -14.802  -1.253   3.199  72.634

Coefficients:
                Estimate Std. Error t value Pr(>|t|)
(Intercept)    5263.1585    15.0667 349.324 4.16e-16 ***
data$Anno        14.3638     8.1282   1.767   0.1205
I(data$Anno^2)  -10.5886     1.1293  -9.376 3.27e-05 ***
I(data$Anno^3)    0.8401     0.4209   1.996   0.0861 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 33.08 on 7 degrees of freedom
Multiple R-squared: 0.9622, Adjusted R-squared: 0.946
F-statistic: 59.44 on 3 and 7 DF, p-value: 2.403e-05


E l'equazione sarebbe:



In quest'ultimo caso però i coefficienti beta1 e beta3 non sono significativi, quindi il modello migliore è il polinomio di 2° grado.

Il problema più grande è adesso quello di rappresentare graficamente quanto ottenuto. In R difatti non esiste (per quanto ne so io) una funzione per plottare i polinomi trovati fittando i dati. Dobbiamo quindi procedere con artefatti grafici comunque validi, ma alquanto laboriosi.

Innanzitutto mettiamo in grafico i valori, con il comando visto prima. Questa volta visualizzeremo solamente le linee e non i punti, per comodità grafica:

plot(data$Anno, data$Popolazione, type="l", lwd=3)

Ora aggiungiamo a questo grafico l'andamento del polinomio di 2° grado in questo modo:

points(data$Anno, predict(fit2), type="l", col="red", lwd=2)

In pratica la funzione predict() non fa altro che sostituire tutti i valori contenuti nella colonna data$Anno nel polinomio di 2° grado che abbiamo calcolato. Quindi le varie coppie di valori vengono unite da segmenti. Quello che otteniamo (vedi figura sotto) è quindi l'unione dei punti virtuali, ossia del valori y = f(x). Non è disegnato il grafico nel continuo, ma nel discreto. Con pochi valori, questo metodo è fortemente debilitante.



Proviamo ad aggiungere anche il grafico del polinomio di 3° grado:

points(data$Anno, predict(fit3), type="l", col="blue", lwd=2)



Come vedete i due modelli hanno andamenti molto simili.

Se volessimo invece ottenere il grafico continuo delle funzioni ottenute, dobbiamo procedere in questo modo. Innanzitutto si crea la funzione polinomiale trovata:

andamento2 <- function(x) fit2$coefficient[3]*x^2 + fit2$coefficient[2]*x + fit2$coefficient[1]

La funzione function(x) precede la formula che è quella del polinomio di 2° grado. fit2$coefficient[3] è beta2, il coefficient[2] è beta1 e coefficient[1] è beta0.

A questo punto mettiamo in grafico solo i punti, e poi visualizziamo la funzione creata con la sintassi curve(x):

plot(data$Anno, data$Popolazione, type="p", lwd=3)
andamento2 <- function(x) fit2$coefficient[3]*x^2 + fit2$coefficient[2]*x + fit2$coefficient[1]
curve(andamento2, col="red", lwd=2)


I punti però scompaiono, ma possiamo riposizionarli con il comando points:

points(data$Anno, data$Popolazione, type="p", lwd=3)

Una nota: è necessario seguire l'ordine dei comandi come ve l'ho descritto, altrimenti la funzione curve genera un grafico errato. Quindi riassumendo i comandi per ottenere la funzione continua, e i punti sperimentali sullo stesso grafico sono i seguenti:

plot(data$Anno, data$Popolazione, type="p", lwd=3)
andamento2 <- function(x) fit2$coefficient[3]*x^2 + fit2$coefficient[2]*x + fit2$coefficient[1]
curve(andamento2, col="red", lwd=2)
points(data$Anno, data$Popolazione, type="p", lwd=3)


Il grafico che otteniamo è il seguente:



Adesso cerchiamo il grafico del polinomio di 3° grado:

plot(data$Anno, data$Popolazione, type="p", lwd=3)
andamento3 <- function(x) fit3$coefficient[4]*x^3 + fit3$coefficient[3]*x^2 + fit3$coefficient[2]*x + fit3$coefficient[1]
curve(andamento3, col="red", lwd=2)
points(data$Anno, data$Popolazione, type="p", lwd=3)




Ribadisco quindi che per ottenere un output corretto è necessario seguire la serie di istruzioni viste finora nell'ordine con cui vi è stato presentato, altrimenti i grafici risulteranno alterati.
Nonostante queste piccole difficoltà, il risultato è pur sempre notevole.

E adesso un piccolo esercizio per voi :-D
Con la serie di dati che qui vi presento, scrivete nei commenti il grado del polinomio e il polinomio stesso che ho trovato per ottenere la curva rossa (dovrete andare a tentativi, così potrete esercitarvi con le istruzioni viste finora, fino a generare un grafico simile al mio!).

x <- c(1:18)
y <- c(7, 6, 4, 3, 4, 5, 4, 3, 2, 3, 4, 5, 6, 4, 2, 3, 5, 7)


mercoledì 22 luglio 2009

Analisi del trend dei visitatori di un sito con i dati di Google Analytics

Questo post è un riassunto di due miei post precedenti sull'analisi del trend con il test di Cox-Stuart e sulla regressione lineare semplice. L'obiettivo che ci proponiamo è quello di valutare il trend del numero di visite ricevute da un sito in un ampio arco temporale. Personalmente utilizzo Google Analytics, perchè questo strumento ci permette di salvare i vari reportage in file Excel. Vediamo punto per punto come salvare il reportage, quindi come importare i dati da Excel in R, ed infine come stimare se il numero di visitatori giornalieri segue un trend in aumento o in diminuzione.

Cominciamo col creare un rapporto ad hoc in Google Analytics. Dopo aver eseguito il login, selezioniamo l'intervallo di date che vogliamo analizzare. Quindi clickiamo su Visitatori.



Nel sottomenù che si apre clickiamo su Visite. A questo punto possiamo salvare il report, clickando su Esporta e quindi clickando su CSV per Excel (è importante salvare in questo formato, e non semplicemente in CSV, altrimenti importare i dati da Excel ad R risulterà difficoltoso).



Salviamo il file CSV, e apriamolo con Excel. Ecco come ci appare (ho preso solo alcune celle):



A questo punto importiamo i dati in R. Con mia grande gioia, ho scoperto che importare i dati da Excel in R è estremamente semplice. E' sufficiente selezionare la colonna (o le colonne) di nostro interesse (nel nostro caso la colonna Visite) e copiare negli appunti (ricordiamoci di selezionare la cella Visite, perchè ci tornerà utile):



Fatto ciò, apriamo R e incolliamo i nostri valori in maniera estremamente semplice con un solo comando:

> myvisit <- read.delim("clipboard")

In che formato R interpreta i dati che gli sono stati forniti? Osserviamo il contenuto di myvisit (qui ci sono solo i primi 10 valori, per comodità):

> myvisit
    Visite
1        6
2       10
3        7
4        8
5        7
6        6
7        4
8        9
9        8
10      11


E' un dataframe a una sola colonna, il cui nome è Visite (per questo ho specificato di copiare anche la cella Visite da Excel).

A questo punto abbiamo in R e possiamo procedere con l'analisi del trend, nei due modi proposti: tramite un test di Cox-Stuart e tramite l'analisi della retta di regressione.

La funzione per effettuare il test di Cox-Stuart è disponibile in questo post. Qui mi limiterò esclusivamente ad utilizzarla. Innanzitutto dobbiamo convertire il dataframe in un formato che può essere letto dalla funzione cox.stuart.test, in questo modo:

> visite <- c(myvisit$Visite)

Ho creato in questo modo un vettore (visite) che contiene tutti i dati ordinati che si trovavano nella colonna Visite del dataframe myvisit. Adesso effettuiamo un test di Cox-Stuart:

> cox.stuart.test(visite)

        Cox-Stuart test for trend analysis

data:
Trend in aumento - p-value = 0


L'output è molto esplicativo: è stato rilevato un trend in aumento delle visite, altamente significativo (essendo p-value < 0.5). Osserviamo qui che p-value = 0. In realtà questa è una approssimazione, perchè eseguendo i calcoli singolarmente (come descritto nel post a riguardo), otterremo p-value = 2.375132e-08, cioè un valore estremamente piccolo.
Riassumendo quindi, secondo il test di Cox-Stuart esiste un significativo trend in aumento delle visite.




Se non siamo soddisfatti o sicuri di questo risultato, possiamo prendere in considerazione il coefficiente angolare della retta di regressione.
Innanzitutto può essere comodo visualizzare i risultati. Il vettore visite contiene le visite giornaliere al sito. Ora creiamo un vettore ordinato dei giorni considerati, della stessa lunghezza del vettore visite:

> giorni <- c(1 : length(visite))

Creiamo un plot con i dati:

> plot(giorni, visite, type="b")

Scegliendo type="b" visualizzo punti e linee, come appare in figura:



Da questo plot non è facile osservare un eventuale trend sull'andamento delle visite. Possiamo però fare un'analisi di regressione. Valutando il segno del coefficiente angolare della retta, possiamo stimare se il trend è in aumento o in calo:

> fit <- lm(visite ~ giorni)
> summary(fit)

Call:
lm(formula = visite ~ giorni)

Residuals:
    Min      1Q  Median      3Q     Max
-33.254  -7.628  -1.253   7.497 106.559

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) 25.93965    1.62693  15.944  < 2e-16 ***
giorni       0.06251    0.01015   6.161 2.56e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 13.5 on 275 degrees of freedom
Multiple R-squared: 0.1213, Adjusted R-squared: 0.1181
F-statistic: 37.96 on 1 and 275 DF, p-value: 2.559e-09


Il coefficiente angolare ha valore b = 0.06251. Esso ha quindi segno positivo (pur essendo piccolo in valore assoluto), quindi è lecito pensare ad un trend in aumento. Il valore della statistica t-test sul coefficiente angolare, e il suo relativo p-value indicano entrambi che esso è significativo. Possiamo quindi affermare che c'è un trend in aumento.

Possiamo infine visualizzare la retta di regressione direttamente sul grafico plot ottenuto in precedenza in questo modo:

> plot(giorni, visite, type="b")
> abline(fit, col="red", lwd=3)


Oppure analogamente in questa maniera:

> plot(giorni, visite, type="b")
> abline(lm(visite ~ giorni), col="red", lwd=3)


Il comando abline permette di aggiungere una retta definita dall'equazione contenuta, direttamente sul grafico che stavamo visualizzando; il parametro col specifica il colore, e il parametro lwd specifica lo spessore della linea. Osserviamo ora il grafico:



E' evidente come ci sia un trend in aumento, sia pur relativamente basso.

martedì 21 luglio 2009

Regressione logistica semplice su variabili qualitative dicotomiche

In questo post vedremo brevemente come realizzare un modello di regressione logistica qualora si disponga di variabili categoriali, o qualitative, organizzate in tabelle di contingenza a doppia entrata. In questo modello la variabile dipendente Y è una variabile bernoulliana, che può assumere i valori 0 oppure 1. La probabilità che assuma valore 1 dipende dai regressori, ossia dalle variabili indipendenti che si prendono in considerazione; anch'esse sono variabili dicotomiche.

In generale, la probabilità che sia Y = 1 in funzione dei regressori è:



Il nostro obiettivo è quello di stimare i valori dei parametri beta.

Cominciamo ad esaminare un modello di regressione logistica semplice (un solo regressore). Sia la variabile dipendente che quella indipendente può assumere valori 0 o 1.

Si consideri il seguente esempio. La tabella sottostante riporta i risultati di uno studio sul reflusso gastroesofageo. Si vuole valutare quanto la presenza di una fattore di stress possa influenzare l'insorgenza di questa patologia.



Innanzitutto importiamo i valori in R. Dobbiamo creare una tabella a doppia entrata; procediamo in questo modo:


reflusso <- matrix(c(251,131,4,33), nrow=2)
colnames(reflusso) <- c("reflNO", "reflSI")
rownames(reflusso) <- c("stressNO", "stressSI")
table <- as.table(reflusso)


Osserviamo la tabella richiamando la variabile:


table

reflNO reflSI
stressNO 251 4
stressSI 131 33


Ora aggiustiamo i dati per effettuare la regressione logistica. Dobbiamo creare un data frame:


dft <- as.data.frame(table)
dft

Var1 Var2 Freq
1 stressNO reflNO 251
2 stressSI reflNO 131
3 stressNO reflSI 4
4 stressSI reflSI 33


Possiamo ore procedere fittando i dati, e quindi effetuando la regressione logistica:


fit <- glm(Var2 ~ Var1, weights = Freq, data = dft, family = binomial(logit))
summary(fit)

Call:
glm(formula = Var2 ~ Var1, family = binomial(logit), data = dft,
weights = Freq)

Deviance Residuals:
1 2 3 4
-2.817 -7.672 5.765 10.287

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.1392 0.5040 -8.213 < 2e-16 ***
Var1stressSI 2.7605 0.5403 5.109 3.23e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 250.23 on 3 degrees of freedom
Residual deviance: 205.86 on 2 degrees of freedom
AIC: 209.86

Number of Fisher Scoring iterations: 6



Questo che vediamo è l'output della funzione. Anzitutto commentiamo il codice per effettuare la regressione. La regressione logistica viene chiamata imponendo la famiglia: family = binomial(logit). Il codice Var2 ~ Var1 significa che vogliamo creare un modello che ci spieghi la variabile Var2 (presenza o assenza del reflusso) in funzione della variabile Var1 (presenza o assenza di eventi stressogeni). In pratica Var2 è la variabile indipendente Y e Var1 è la variabile dipendente X (il regressore). Fornita la formula da analizzare, si specifica il peso di ciascuna variabile, dati contenuti nella colonna Freq del dataframe dft (quindi si scrive weights = Freq e data = dft per specificare la posizione dove sono contenuti i valori). Con summary(fit) si ottengono una serie di dati utili.

I valori dei parametri $\beta_0$ e $\beta_1$ sono rispettivamente i valori (intercept) e Var1stress1. Possiamo quindi scrivere il nostro modello empirico:



La variabile indipendente x può assumere valore zero oppure uno. Se assume valore 0 (quindi in assenza di eventi stressanti), allora la probabilità di avere il reflusso è pari a:



Se sono presenti eventi stressanti (x=1), la probabilità di avere reflusso è:



Le formula per gli odds sono:



Possiamo infine calcolare l'odd ratio OR:



Una persona che ha vissuto un evento stressante ha una propensione a sviluppare il reflusso gastroesofageo 15.807 volte più grande della persona che non ha subito eventi stressanti.

Gli odds, le probabilità, e l'OR, possono essere rapidamente calcolati in R ricordando che:

fit$coefficient[1] = (intercept)
fit$coefficient[2] =

Inoltre:

summary(fit)$coefficient[1,2] = standard error del coefficiente $\beta_0$
summary(fit)$coefficient[2,2] = standard error del coefficiente $\beta_1$

E quindi abbiamo:


pi0 <- exp(fit$coefficient[1]) / (1 + exp(fit$coefficient[1]))
pi1 <- exp(fit$coefficient[1] + fit$coefficient[2]) / (1 + exp(fit$coefficient[1]+fit$coefficient[2]))

odds0 <- pi0 / (1 - pi0)
odds1 <- pi1 / (1 - pi1)

OR <- odds1 / odds0

#stesso risultato con:
OR <- exp(fit$coefficient[2])

#l'intervallo di confidenza per l'odd-ratio è:
ORmin <- exp( fit$coefficient[2] - qnorm(.975) * summary(fit)$coefficient[2,2] )

ORmax <- exp( fit$coefficient[2] + qnorm(.975) * summary(fit)$coefficient[2,2] )


La più nota formula per l'odd-ratio è la seguente (stesso risultato):



che in R calcoliamo così:


OR <- (table[1,1]*table[2,2]) / (table[1,2]*table[2,1])


L'acronimo AIC sta per Akaike's information criterion. Questo parametro non ci fornisce alcun dato sul modello appena creato. Esso può però essere utile nel confrontare questo modello con altri eventualmente presi in considerazione (il modello con AIC più basso è il migliore).

La regressione logistica multipla viene trattata in questo post

sabato 18 luglio 2009

Analisi del trend con il test di Cox-Stuart in R

Il Cox-Stuart test è definito un test poco potente (potenza pari a 0.78), ma molto robusto per l'analisi del trend (trend analysis). E' quindi applicabile ad una grande varietà di situazioni, per avere un'idea dell'andamento dei valori ottenuti. Il metodo proposto è fondato sulla distribuzione binomiale. In R non esiste una funzione per eseguire un test di Cox-Stuart, pertanto vediamo i passaggi logici che sono alla base del test, e infine potremo scrivere noi stessi la funzione adatta.

Si vuole valutare se esiste un trend in aumento o in diminuzione per una serie di rilevazioni effettuate misurando il numero di clienti al giorno in un ristorante nell'arco di 15 giorni. Questi sono i valori:
Clienti: 5, 9, 12, 18, 17, 16, 19, 20, 4, 3, 18, 16, 17, 15, 14


Per eseguire il test di Cox-Stuart, il numero di rilevazioni deve essere pari. Nel nostro caso abbiamo 15 rilevazioni. Eliminiamo pertanto la rilevazione in posizione (N+1)/2 (nel nostro caso la misura con valore 20):

> clienti = c(5, 9, 12, 18, 17, 16, 19, 20, 4, 3, 18, 16, 17, 15, 14)
>
> length(clienti)
[1] 15
> clientipari = clienti[ -(length(clienti)+1)/2]
> length(clientipari)
[1] 14


Ora abbiamo 14 misure e possiamo quindi procedere. Dividiamo le misure in due vettori, contenenti il primo la prima metà delle misure, e il secondo la seconda metà:

> meta1 = clientipari[1:7]
> meta2 = clientipari[8:14]
> meta1
[1] 5 9 12 18 17 16 19
> meta2
[1] 4 3 18 16 17 15 14


Adesso sottraiamo valore per valore il contenuto dei due vettori:

> differenze = meta1 - meta2
> differenze
[1] 1 6 -6 2 0 1 5


Da questo vettore, estraiamo i segni di ogni valore:

> segni = sign(differenze)
> segni
[1] 1 1 -1 1 0 1 1


Osserviamo che una differenza ha dato valore 0 e pertanto anche nel vettore con i segni vi è un valore pari a 0. Questo deve essere eliminato:

> segni = segni[ segni != 0 ]
> segni
[1] 1 1 -1 1 1 1


Abbiamo ottenuto 6 differenze, e quindi 6 segni. A questo punto dobbiamo contare quanti sono i segni positivi e quanti i segni negativi. In questo caso, si può facilmente contare a mano. Se le rilevazioni sono però numerose possiamo procedere così:

> pos = segni[segni > 0]
> neg = segni[segni < 0]
> length(pos)
[1] 5
> length(neg)
[1] 1


A questo punto si sceglie il numero di segni minore. In questo caso i segni negativi, pari a 1. Si calcola la probabilità (secondo una distribuzione binomiale) di avere x=1 successo su N=6 rilevazioni, con probabilità al singolo gioco pari a 0.5:

> pbinom(1, 6, 0.5)
[1] 0.109375


Il valore così calcolato è superiore a 0.05 (significatività del 95%). Pertanto non esiste un trend significativo (che sarebbe stato in diminuzione, essendo il numero di segni negativi minore).
Se il valore fosse stato inferiore a 0.05, avremmo accettato come significativo il trend.

Ecco qui di seguito il codice per effettuare immediatamente il test di Cox-Stuart, scritto da me.

cox.stuart.test =
function (x)
{
method = "Cox-Stuart test for trend analysis"
leng = length(x)
apross = round(leng) %% 2
if (apross == 1) {
elimina = (length(x)+1)/2
x = x[ -elimina ]
}
meta = length(x)/2
x1 = x[1:meta]
x2 = x[(meta+1):(length(x))]
difference = x1-x2
segni = sign(difference)
segnitolti = segni[segni != 0]
pos = segni[segni>0]
neg = segni[segni<0]
if (length(pos) < length(neg)) {
prop = pbinom(length(pos), length(segnitolti), 0.5)
names(prop) = "Trend in aumento - p-value"
rval <- list(method = method, statistic = prop)
class(rval) = "htest"
return(rval)
}
else {
prop = pbinom(length(neg), length(segnitolti), 0.5)
names(prop) = "Trend in diminuzione - p-value"
rval <- list(method = method, statistic = prop)
class(rval) = "htest"
return(rval)
}
}


Scritta la funzione in R, vediamo subito il suo output inserendo gli stessi dati analizzati poc'anzi:

> clienti = c(5, 9, 12, 18, 17, 16, 19, 20, 4, 3, 18, 16, 17, 15, 14)
> cox.stuart.test(clienti)

Cox-Stuart test for trend analysis

data:
Probabilità in diminuzione = 0.1094


L'output è il valore della probabilità, del quale bisogna valutare la significatività confrontandolo con il valore 0.05 (0.01, o 0.005, a seconda delle necessità), come visto prima.

venerdì 17 luglio 2009

ANOVA a due vie

L'analisi della varianza ad una via è utile per verificare contemporaneamente se le medie di più gruppi sono uguali. Ma questa analisi può risultare poco utile, ai fini di problemi più complessi. Ad esempio può essere necessario prendere in considerazione due fattori di variabilità, per verificare se le medie tra i gruppi dipendono dal gruppo di classificazione ("zone") o dalla seconda variabile che si va a considerare ("blocco"). In questo caso si ricorre ad una analisi della varianza a due vie (ANOVA a due vie, two-way ANOVA).

Cominciamo subito con un esempio, così da rendere più facile la comprensione di questo metodo statistico. I dati raccolti vengono organizzati in tabelle a doppia entrata.

Il direttore di una società ha raccolto le entrate (in migliaia di dollari) per 5 anni e in base al mese. Si vuole verificare se le entrate dipendono dall'annata e/o dal mese, oppure se sono indipendenti da questi due fattori.

Concettualmente il problema potrebbe essere risolto con un'ANOVA orizzontale, ed un ANOVA verticale. Si verifica cioè se le entrate medie classificate per anno siano uguali, e se sono uguali le entrate medie classificate per mese. Ciò richiederebbe numerosi calcoli, e pertanto si ricorre all'analisi ANOVA a due vie che offre il risultato istantaneamente.
Questa è la tabella delle entrate classificate per anno e per mese:



Come per l'ANOVA a una via, anche qui l'obiettivo è quello di strutturare un test F di Fisher per valutare la significatività della variabile mesi e della variabile anno, e stabilire se in base a uno (o entrambi) di questi criteri di classificazione dipendono le entrate.
In R la sintassi per l'ANOVA a due vie è la seguente. Innanzitutto si crea un array contenente tutti i valori tabulati, trascritti per riga:

> entrate = c(15,18,22,23,24, 22,25,15,15,14, 18,22,15,19,21,
+ 23,15,14,17,18, 23,15,26,18,14, 12,15,11,10,8, 26,12,23,15,18,
+ 19,17,15,20,10, 15,14,18,19,20, 14,18,10,12,23, 14,22,19,17,11,
+ 21,23,11,18,14)
>


A questo punto si classificano i valori così inseriti. In base ai mesi si crea un fattore a 12 livelli (numero righe) con 5 ripetizioni (numero colonne) in questa maniera:
> mesi = gl(12,5)

In base agli anni si crea un fattore a 5 livelli (numero colonne) e 1 ripetizione; si esplicita il numero di variabili inserite (la lunghezza dell'array entrate):
> anni = gl(5, 1, 60)

A questo punto è possibile fittare il modello lineare e produrre la tabella ANOVA:

> fit = aov(entrate ~ mesi + anni)
>
> anova(fit)

Analysis of Variance Table

Response: entrate
          Df Sum Sq Mean Sq F value Pr(>F)
mesi      11 308.45   28.04  1.4998 0.1660
anni       4  44.17   11.04  0.5906 0.6712
Residuals 44 822.63   18.70


Interpretiamo adesso i risultati.
La significatività della differenza tra mesi è: F = 1.4998. Questo valore è inferiore al valore tabulato e difatti p-value > 0.05. Quindi si accetta l'ipotesi nulla<: le medie valutate in base ai mesi sono uguali; quindi la variabile "mesi" non influisce sulle entrate.

La significatività della differenza tra anni è: F = 0.5906. Questo valore è inferiore al valore tabulato e difatti p-value > 0.05. Quindi si accetta l'ipotesi nulla: le medie valutate in base agli anni sono uguali; quindi la variabile "anni" non influisce sulle entrate.

martedì 7 luglio 2009

Regressione lineare

Si ricorre alla analisi della regressione quando dai dati campionari si vuole ricavare un modello statistico che predica i valori di una variabile (Y) detta dipendente a partire dai valori di un'altra variabile (X) detta indipendente.
La regressione lineare, che rappresenta la relazione più semplice e frequente tra due variabili quantitative, può essere positiva (all'aumento dei valori di una variabile corrisponde un aumento anche nell'altra) o negativa (all'aumento dell'una corrisponde una diminuzione dell'altra): tale relazione è indicata dal segno del coefficiente b.
Ricordo le formule per calcolare il coefficiente angolare e l'intercetta:



Per costruire la retta che descrive la distribuzione dei punti, ci si può riferire a diversi principi. Il più comune è il metodo dei minimi quadrati (least squares, o Model 1), ed è il metodo utilizzato dal software statistico R.

Supponiamo di voler ricavare una relazione lineare tra il peso (kg) e l'altezza (cm) di 10 individui.
Altezza: 175, 168, 170, 171, 169, 165, 165, 160, 180, 186
Peso: 80, 68, 72, 75, 70, 65, 62, 60, 85, 90


Il primo problema che si pone è quello di decidere quale sia la variabile dipendente Y e quale la variabile indipendente X. In generale, la variabile indipendente è quella non affetta da errore durante la misura (o affetta da errore casuale); mentre la variabile dipendente è quella affetta da errore, e di cui si vuole stimare una relazione. Nel nostro caso possiamo assumere che la variabile Peso sia la variabile indipendente (X), e la variabile Altezza quella dipendente (Y).
Quindi il nostro problema è quello di cercare una relazione lineare (in termini spiccioli, una formula) che ci permetta di calcolare l'altezza, essendo noto il peso di un individuo. La formula più semplice è quella di una generica retta del tipo Y = a + bX. In R si calcolano i due parametri procedendo in questo modo:

> altezza = c(175, 168, 170, 171, 169, 165, 165, 160, 180, 186)
> peso = c(80, 68, 72, 75, 70, 65, 62, 60, 85, 90)
>
> modello = lm(formula = altezza ~ peso, x=TRUE, y=TRUE)
> modello

Call:
lm(formula = altezza ~ peso, x = TRUE, y = TRUE)

Coefficients:
(Intercept)         peso
   115.2002       0.7662


Anzitutto, notiamo che per scrivere correttamente la funzione, è importante l'ordine delle variabili. In formula = altezza ~ peso, si è seguita la sintassi generale formula = Y ~ X, e questa sintassi è stata confermata esplicitando x=TRUE e y=TRUE.
L'output della funzione è rappresentato dai due parametri a e b: a=115.2002 (l'intercetta), b=0.7662 (il coefficiente angolare).




Il semplice calcolo della retta non è però sufficiente. Occorre valutare la significatività della retta, ossia se il coefficiente angolare b si discosta da zero in modo significativo. Il test può essere effettuato sia mediante il test F di Fisher, sia con il test t di Student.

Ricordiamo le ipotesi da verificare:



In R entrambi possono essere richiamati molto rapidamente. Ecco come:

> mod = lm(altezza ~ peso)
> summary(mod)

Call:
lm(formula = altezza ~ peso)

Residuals:
    Min      1Q Median      3Q     Max
-1.6622 -0.9683 -0.1622  0.5679  2.2979

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 115.20021    3.48450   33.06 7.64e-10 ***
peso          0.76616    0.04754   16.12 2.21e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.405 on 8 degrees of freedom
Multiple R-squared: 0.9701, Adjusted R-squared: 0.9664
F-statistic: 259.7 on 1 and 8 DF, p-value: 2.206e-07


Come possiamo vedere, l'output è molto ricco. Innanzitutto osserviamo che sono stati forniti i valori dell'intercetta e del coefficiente angolare anche in questo modo (io personalmente preferisco ricavarli col metodo visto prima, per rapidità di interpretazione). Viene inoltre il risultato "Multiple R-squared: 0.835": questo numero descrive la bontà del modello ($R^2$ è il coefficiente di determinazione), ossia quanto il modello trovato spiega i dati campionari. In questo caso il modello è adeguato a descrivere l'83.5% dei dati ricavati.
Il test t di Student relativo al coefficiente angolare ha in questo caso valore 16.12 (quello relativo alla variabile peso); il risultato della statistica F di Fisher è invece 259.7 (è lo stesso valore che si otterrebbe eseguendo una ANOVA sui dati, che possiamo guardare richiamandolo così: anova(mod)). In entrambi i casi questi valori sono ampiamente superiori ai valori tabulati (e difatti i relativi p-value sono di molto inferiori a 0.05), e pertanto viene rifiutata l'ipotesi nulla (b = 0): quindi la regressione è significativa (il valore del coefficiente angolare così calcolato è statisticamente differente da zero).

domenica 5 luglio 2009

Tabella di contingenza per lo studio della correlazione tra variabili qualitative

Qualora si disponga di variabili qualitative, è possibile verificare la correlazione studiando una tabella di contingenza R x C.

Un casinò vuole studiare la correlazione tra la modalità di gioco e il numero di vincitori per classi di età, per verificare se il numero di vincitori dipende dal tipo di gioco che si è scelto di fare, in base all'esperienza. Si dispone dei seguenti dati:
Classi di età: 20-30, 31-40, 41-50
Roulette: 44, 56, 55
Black-jack: 66, 88, 23
Poker: 15, 29, 45


In R dobbiamo anzitutto costruire una matrice con i dati raccolti:

tabella = matrix(c(44,56,55, 66,88,23, 15,29,45), nrow=3, byrow=TRUE)

La sintassi corretta della funzione matrix prevede di esprimere il numero di righe (nrow=3), e di specificare in che ordine sono stati inseriti i dati (in questo caso sono stati ordinati per riga, quindi byrow=TRUE).
Calcoliamo ora il chi-quadro:

> chisq.test(tabella)

Pearson's Chi-squared test

data: tabella
X-squared = 46.0767, df = 4, p-value = 2.374e-09


Il valore X-squared è di molto superiore al valore del chi-quadro-tabulato:

> qchisq(0.950, 4)
[1] 9.487729


Pertanto rifiuto l'ipotesi nulla: esiste correlazione tra i valori analizzati (del resto p-value < 0.05: ci porta a rifiutare l'ipotesi H0, quindi la statistica test è significativa).

Metodi non parametrici per lo studio della correlazione: test di Spearman e di Kendall

Abbiamo visto nel post precedente, in che modo studiare la correlazione tra variabili che seguano una distribuzione gaussiana con il metodo di Pearson. Nel caso in cui non sia possibile assumere che i valori seguano distribuzioni gaussiane, sono a disposizione due metodi non parametrici: il test rho di Spearman e il test tau di Kendall.

Si vuole studiare la produttività di diversi tipi di macchinari e la soddisfazione degli operatori nel loro utilizzo (espressi con un numero da 1 a 10). Questi sono i valori:
Produttività: 5, 7, 9, 9, 8, 6, 4, 8, 7, 7
Soddisfazione: 6, 7, 4, 4, 8, 7, 3, 9, 5, 8


Cominciamo ad utilizzare per primo il test Rho di Spearman:

> a = c(5, 7, 9, 9, 8, 6, 4, 8, 7, 7)
> b = c(6, 7, 4, 4, 8, 7, 3, 9, 5, 8)
>
> cor.test(a, b, method="spearman")

Spearman's rank correlation rho

data: a and b
S = 145.9805, p-value = 0.7512
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
0.1152698


La statistica test ci dà come risultato rho = 0.115, il che indica una bassa correlazione (non parametrica) tra le due serie di valori.
Il p-value > 0.05 ci indica che il test non è statisticamente significativo (ossia il valore di Rho campionario non può essere esteso a tutta la popolazione.

Verifichiamo adesso gli stessi dati con il test Tau di Kendall:

> a = c(5, 7, 9, 9, 8, 6, 4, 8, 7, 7)
> b = c(6, 7, 4, 4, 8, 7, 3, 9, 5, 8)
>
> cor.test(a, b, method="kendall")

Kendall's rank correlation tau

data: a and b
z = 0.5555, p-value = 0.5786
alternative hypothesis: true tau is not equal to 0
sample estimates:
tau
0.146385


Anche con il test di Kendall, la correlazione risulta essere molto bassa (tau = 0.146), e p-value > 0.05.

Metodo parametrico per lo studio della correlazione: test R di Pearson

Supponiamo di voler studiare se esiste una correlazione tra 2 serie di dati. Per fare questo ricorriamo al test R di Pearson, che ci fornisce come risultato quanto è "forte" la correlazione e il valore del test t per studiare la significatività del coefficiente R di Pearson.

Un nuovo test per misurare il QI viene sottoposto a 10 volontari. Si vuole verificare se esiste una correlazione tra i risultati al nuovo test sperimentale, e i risultati al classico test. Questi i valori:
Test-vecchio: 15, 21, 25, 26, 30, 30, 22, 29, 19, 16
Test-nuovo: 55, 56, 89, 67, 84, 89, 99, 62, 83, 88


Nel software R è disponibile un'unica funzione, facilmente richiamabile, che ci fornisce direttamente il valore del coefficiente di Pearson e la statistica test per la verifica della significatività del coefficiente:

> a = c(15, 21, 25, 26, 30, 30, 22, 29, 19, 16)
> b = c(55, 56, 89, 67, 84, 89, 99, 62, 83, 88)
>
> cor.test(a, b)

Pearson's product-moment correlation

data: a and b
t = 0.4772, df = 8, p-value = 0.646
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.5174766 0.7205107
sample estimates:
cor
0.166349


Il valore del coefficiente di Pearson è 0.166: è un valore molto basso, che indica una scarsa correlazione tra le variabili.
Controlliamo ora la significatività di R dal valore della statistica test:

> qt(0.950, 8)
[1] 1.859548


Risulta che t-calcolato < t-tabulato; del resto p-value > 0.05. Quindi il coefficiente R di Pearson non è statisticamente significativo.

Confronto tra più gruppi, metodo non parametrico: test di Kruskal-Wallis

Se vogliamo eseguire il confronto tra più gruppi, ma non è possibile eseguire una ANOVA, è possibile ricorrere al test di Kruskal-Wallis, che può essere applicato quando non si può fare l'assunzione che i gruppi seguano una distribuzione normale.
Questo test è del tutto analogo ai test di Wilcoxon per 2 campioni.
Può essere utile verificare l'omogeneità delle varianze di k gruppi di campioni che non seguono una distribuzione normale, applicando il test di Levene non parametrico.

Supponiamo di volere verificare se i seguenti 4 gruppi di valori siano statisticamente simili:
Gruppo A: 1, 5, 8, 17, 16
Gruppo B: 2, 16, 5, 7, 4
Gruppo C: 1, 1, 3, 7, 9
Gruppo D: 2, 15, 2, 9, 7


Per richiamare il test di Kruskal-Wallis è sufficiente inserire i dati, e quindi organizzarli in una lista:

> a = c(1, 5, 8, 17, 16)
> b = c(2, 16, 5, 7, 4)
> c = c(1, 1, 3, 7, 9)
> d = c(2, 15, 2, 9, 7)
>
> dati = list(g1=a, g2=b, g3=c, g4=d)


Se richiamiamo la variabile dati, possiamo vedere come i dati sono stati organizzati.
Ora richiamiamo la funzione:

> kruskal.test(dati)

Kruskal-Wallis rank sum test

data: dati
Kruskal-Wallis chi-squared = 1.9217, df = 3, p-value = 0.5888


Il valore della statistica test è 1.9217. Questo valore prevede già la correzione quando ci sono dei ties (ripetizioni). Il p-value è maggiore di 0.05; inoltre il valore della statistica test è inferiore al chi-quadro-tabulato:

> qchisq(0.950, 3)
[1] 7.814728


La conclusione è quindi che accetto l'ipotesi H0: le medie dei 4 gruppi sono statisticamente uguali.

sabato 4 luglio 2009

Confronto tra più gruppi, metodo parametrico: analisi della varianza ANOVA

Analisi della varianza: ANOVA

Il metodo ANOVA permette di effettuare il confronto tra più gruppi, con metodo parametrico (supponendo cioè che i vari gruppi seguano una distribuzione gaussiana).
Procediamo con il seguente esempio:

Il dirigente di una catena di supermercati vuole verificare se il consumo in KiloWatt di 4 suoi negozi siano tra loro uguali. Raccoglie i dati alla fine di ogni mese, per 6 mesi. I risultati sono i seguenti:
Supermercato A: 65, 48, 66, 75, 70, 55
Supermercato B: 64, 44, 70, 70, 68, 59
Supermercato C: 60, 50, 65, 69, 69, 57
Supermercato D: 62, 46, 68, 72, 67, 56


Per procedere con la verifica ANOVA, occorre dapprima verificare l'omoschedasticità (ossia effettuare test per l'omogeneità delle varianze). Il software R mette a disposizione due test nel pacchetto standard: il test di Bartlett e il test di Fligner-Killeen; numerosi altri test per testare l'omogeneità delle varianze (test di Hartley, di Cochran, di Levene, eccetera) vengono discussi in questo post.




Cominciamo con il test di Bartlett.

Per prima cosa creiamo 4 array con i valori:


a = c(65, 48, 66, 75, 70, 55)
b = c(64, 44, 70, 70, 68, 59)
c = c(60, 50, 65, 69, 69, 57)
d = c(62, 46, 68, 72, 67, 56)


Adesso concateniamo i 4 valori in un unico vettore:


dati = c(a, b, c, d)


Ora su questo unico contenitore in cui sono salvati tutti i dati, creiamo i 4 livelli:


gruppi = factor(rep(letters[1:4], each = 6))


Possiamo osservare il contenuto del vettore gruppi semplicemente digitando gruppi + [invio].

A questo punto avviamo il test di Bartlett:


bartlett.test(dati, gruppi)

Bartlett test of homogeneity of variances

data: dati and gruppi
Bartlett's K-squared = 0.4822, df = 3, p-value = 0.9228


La funzione ci ha fornito il valore della statistica test, nonché il p-value. Possiamo affermare che le varianze sono omogenee essendo p-value > 0.05. In alternativa possiamo confrontare il Bartlett's K-squared con il valore del chi-quadro-tabulato; calcoliamo quest'ultimo valore, assegnando il valore di alpha e i gradi di libertà:


qchisq(0.950, 3)
[1] 7.814728


Come vedete, chi-quadro-tabulato è maggiore di Bartlett's K-squared, quindi l'ipotesi di omoschedasticità è vera.




Proviamo a verificare l'omoschedasticità con il test di Fligner-Killeen.
La sintassi è del tutto analoga, quindi procediamo velocemente.


a = c(65, 48, 66, 75, 70, 55)
b = c(64, 44, 70, 70, 68, 59)
c = c(60, 50, 65, 69, 69, 57)
d = c(62, 46, 68, 72, 67, 56)

dati = c(a, b, c, d)

gruppi = factor(rep(letters[1:4], each = 6))

fligner.test(dati, gruppi)

Fligner-Killeen test of homogeneity of variances

data: dati and gruppi
Fligner-Killeen:med chi-squared = 0.1316, df = 3, p-value = 0.9878


Le conclusioni sono analoghe a quelle viste per il test di Bartlett.




Verificata l'omoschedasticità dei 4 gruppi, possiamo procedere con il metodo ANOVA.

Anzitutto organizziamo i valori:


modello = lm(formula = dati ~ gruppi)


Quindi analizziamo col metodo ANOVA:


anova (modello)

Analysis of Variance Table

Response: dati
Df Sum Sq Mean Sq F value Pr(>F)
gruppi 3 8.46 2.82 0.0327 0.9918
Residuals 20 1726.50 86.33


L'output della funzione è la classica tabella ANOVA con i seguenti dati:
Df = gradi di libertà
Sum Sq = devianza (entro gruppi, e residua)
Mean Sq = varianza (entro gruppi, e residua)
F value = è il valore della statistica test, calcolato come (varianza entro gruppi) / (varianza residua)
Pr(>F) = è il p-value

Essendo p-value > 0.05, accettiamo l'ipotesi H0: le medie dei valori sono statisticamente uguali. Per sicurezza possiamo confrontare F-value con il valore F-tabulato (alpha, gradi di libertà del numeratore, gradi di libertà del denominatore):


qf(0.950, 20, 3)
[1] 8.66019


Poiché F-tabulato > F-value-calcolato, posso concludere accettando l'ipotesi H0: le medie sono statisticamente uguali.

Se l'ANOVA fosse risultata significativa (ossia se p-value < 0.05, allora vuol dire che non tutte le medie dei k gruppi considerati sono tra loro uguali. Occorre andare a ricercare quale o quali coppie di gruppi rendono significativa l'ANOVA, ossia si deve cercare quali coppie di medie sono significativamente differenti. Per far questo sono a disposizione gli ANOVA post hoc-tests, di cui si parla in questo post.