R Language
Realización de una prueba de permutación

Una función bastante general

Utilizaremos el conjunto de datos de crecimiento dental incorporado. Nos interesa saber si existe una diferencia estadísticamente significativa en el crecimiento de los dientes cuando a los cobayas se les da vitamina C en comparación con el jugo de naranja.

Aquí está el ejemplo 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

Después de leer el CSV, definimos la función.

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
}

Esta función toma dos vectores, y baraja sus contenidos juntos, luego realiza la función testStat en los vectores barajados. El resultado de teststat se agrega a los trials , que es el valor de retorno.

Lo hace N = 10^5 veces. Tenga en cuenta que el valor N podría muy bien haber sido un parámetro para la función.

Esto nos deja con un nuevo conjunto de datos, trials , el conjunto de medios que podría resultar si realmente no hay relación entre las dos variables.

Ahora para definir nuestra estadística de prueba:

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

Realice la prueba:

result = permutationTest(vec1, vec2, subtractMeans)

Calcule nuestra diferencia de medias observada real:

observedMeanDifference = subtractMeans(vec1, vec2)

Veamos cómo se ve nuestra observación en un histograma de nuestra estadística de prueba.

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

No parece que nuestro resultado observado sea muy probable que ocurra por casualidad ...

Queremos calcular el valor p, la probabilidad del resultado original observado si no hay una relación entre las dos variables.

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

Vamos a desglosar un poco:

result >= (observedMeanDifference)

Creará un vector booleano, como:

FALSE TRUE FALSE FALSE TRUE FALSE ...

Con TRUE cada vez que el valor del result es mayor o igual al valor observedMean .

La mean función interpretará este vector como 1 para TRUE y 0 para FALSE , y nos dará el porcentaje de 1 en la mezcla, es decir, el número de veces que nuestra diferencia de medias del vector barajado supera o iguala lo que observamos.

Finalmente, multiplicamos por 2 porque la distribución de nuestro estadístico de prueba es altamente simétrica, y realmente queremos saber qué resultados son "más extremos" que nuestro resultado observado.

Todo lo que queda es dar salida al valor p, que resulta ser 0.06093939 . La interpretación de este valor es subjetiva, pero diría que parece que la vitamina C promueve el crecimiento de los dientes mucho más que el jugo de naranja.