R Language
Esecuzione di un test di permutazione

Una funzione abbastanza generale

Useremo il set di dati di crescita dei denti integrato. Ci interessa sapere se c'è una differenza statisticamente significativa nella crescita dei denti quando alle cavie viene somministrata vitamina C rispetto al succo d'arancia.

Ecco l'esempio completo:

teethVC = ToothGrowth[ToothGrowth$supp == 'VC',]
teethOJ = ToothGrowth[ToothGrowth$supp == 'OJ',]

permutationTest = function(vectorA, vectorB, testStat){
  N = 10^5
  fullSet = c(vectorA, vectorB)
  lengthA = length(vectorA)
  lengthB = length(vectorB)
  trials <- replicate(N, 
                      {index <- sample(lengthB + lengthA, size = lengthA, replace = FALSE)
                      testStat((fullSet[index]), fullSet[-index])  } )
  trials
}
vec1 =teethVC$len;
vec2 =teethOJ$len;
subtractMeans = function(a, b){ return (mean(a) - mean(b))}
result = permutationTest(vec1, vec2, subtractMeans)
observedMeanDifference = subtractMeans(vec1, vec2)
result = c(result, observedMeanDifference)
hist(result)
abline(v=observedMeanDifference, col = "blue")
pValue = 2*mean(result <= (observedMeanDifference))
pValue

Dopo aver letto nel CSV, definiamo la funzione

permutationTest = function(vectorA, vectorB, testStat){
  N = 10^5
  fullSet = c(vectorA, vectorB)
  lengthA = length(vectorA)
  lengthB = length(vectorB)
  trials <- replicate(N, 
                      {index <- sample(lengthB + lengthA, size = lengthA, replace = FALSE)
                      testStat((fullSet[index]), fullSet[-index])  } )
  trials
}

Questa funzione prende due vettori e mischia i loro contenuti insieme, quindi esegue la funzione testStat sui vettori mescolati. Il risultato di teststat viene aggiunto alle trials , che è il valore di ritorno.

Lo fa N = 10^5 volte. Si noti che il valore N potrebbe essere stato un parametro della funzione.

Questo ci lascia con una nuova serie di dati, trials , l'insieme di mezzi che potrebbe risultare se non ci fosse veramente alcuna relazione tra le due variabili.

Ora per definire la nostra statistica di test:

subtractMeans = function(a, b){ return (mean(a) - mean(b))}

Esegui il test:

result = permutationTest(vec1, vec2, subtractMeans)

Calcola la nostra differenza media osservata effettiva:

observedMeanDifference = subtractMeans(vec1, vec2)

Vediamo come appare la nostra osservazione su un istogramma della nostra statistica di test.

hist(result)
abline(v=observedMeanDifference, col = "blue")

Non sembra che il nostro risultato osservato sia molto probabile che si verifichi per caso ...

Vogliamo calcolare il p-value, il likeliehood del risultato originale osservato se la loro non è una relazione tra le due variabili.

pValue = 2*mean(result >= (observedMeanDifference))

Riassumiamolo un po ':

result >= (observedMeanDifference)

Creerà un vettore booleano, come:

FALSE TRUE FALSE FALSE TRUE FALSE ...

Con TRUE ogni volta il valore del result è maggiore o uguale al valore observedMean .

La mean funzione interpreterà questo vettore come 1 per TRUE e 0 per FALSE e ci darà la percentuale di 1 nel mix, ovvero il numero di volte in cui la differenza media del vettore mescolato ha superato o eguagliato ciò che abbiamo osservato.

Infine, moltiplichiamo per 2 perché la distribuzione della nostra statistica di test è altamente simmetrica e vogliamo veramente sapere quali risultati sono "più estremi" dei nostri risultati osservati.

Tutto ciò che rimane è di generare il valore p, che risulta essere 0.06093939 . L'interpretazione di questo valore è soggettiva, ma direi che sembra che la vitamina C promuova la crescita dei denti molto più di quanto faccia Orange Juice.