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.