Análisis de la Varianza de Friedman con R

Francisco García García (2014-12-12)


1. Objetivo del estudio y diseño del experimento

2. Organización de los datos en R

3. Análisis exploratorio de los datos

4. Detección de diferencias significativas: ANOVA de Friedman

5. Bibliografía y enlaces de interés





1. Objetivo del estudio y diseño del experimento

El objetivo del estudio es evaluar las diferencias del crecimiento celular entre dos grupos experimentales: Wild Type (WT) y un mutante (KO).

Diseño del experimento. Se realizaron un conjunto de experimentos con la misma estructura: 4 puntos de tiempo (24h, 48h, 72h y 96h) donde se cuantificaron el número de células en cada uno de los dos grupos descritos (WT, KO)

2. Organización de los datos en R

El experimento NE77 se excluyó del estudio por presentar diferentes puntos de tiempos

Leemos los datos:

datos <- read.csv("ful.csv",sep = "\t", header = T)

Estructura de los datos:

datos
#>    experiment group    t24     t48     t72     t96
#> 1        NE61    WT 110000  250000  355000  870000
#> 2        NE61    KO 115000  185000  446666 1085000
#> 3        NE62    WT 155000  285000  260000  580000
#> 4        NE62    KO 125000  245000  290000  600000
#> 5        NE63    WT 185000 1540000 2445000 2020000
#> 6        NE63    KO 190000  610000 1990000 2310000
#> 7        NE63    KO 110000  805000 2680000 1500000
#> 8        NE64    WT 290000  185000  420000  480000
#> 9        NE64    KO 170000  235000  410000  473333
#> 10       NE78    WT   9250   16000   20750   36000
#> 11       NE78    KO  27500   33500   40000   75500
#> 12       NE79    WT   7000   45000   56500   82000
#> 13       NE79    KO   7500   76750  123500  291000
#> 14       NE80    WT  18500   60500  154750  200000
#> 15       NE80    WT   9250   34750  133000  190000
#> 16       NE80    KO   4500   16750   90500  118000
#> 17       NE80    KO   5500   18250  102000  140000

Dimensión de la matriz de datos (número de filas x número de columnas):

dim(datos)
#> [1] 17  6



3. Análisis exploratorio de los datos

3.1. Exploramos la muestra determinando algunos estadísticos descriptivos:

summary(datos)
#>  experiment group       t24              t48               t72         
#>  NE61:2     KO:9   Min.   :  4500   Min.   :  16000   Min.   :  20750  
#>  NE62:2     WT:8   1st Qu.:  9250   1st Qu.:  34750   1st Qu.: 102000  
#>  NE63:3            Median :110000   Median : 185000   Median : 260000  
#>  NE64:2            Mean   : 90529   Mean   : 273029   Mean   : 589274  
#>  NE78:2            3rd Qu.:155000   3rd Qu.: 250000   3rd Qu.: 420000  
#>  NE79:2            Max.   :290000   Max.   :1540000   Max.   :2680000  
#>  NE80:4                                                                
#>       t96         
#>  Min.   :  36000  
#>  1st Qu.: 140000  
#>  Median : 473333  
#>  Mean   : 650049  
#>  3rd Qu.: 870000  
#>  Max.   :2310000  
#> 

¿Hay mucha variabilidad entre los 16 experimentos? Sí. Lo comprobamos determianando las desviaciones típicas en cada punto de tiempo y grupo experimental:

desviaciones_wt <- apply(datos[datos$group=="WT",c("t24", "t48", "t72", "t96")], 2, sd)
desviaciones_wt
#>      t24      t48      t72      t96 
#> 105681.7 510869.8 805809.8 655089.9
desviaciones_ko <- apply(datos[datos$group=="KO",c("t24", "t48", "t72", "t96")], 2, sd)
desviaciones_ko
#>       t24       t48       t72       t96 
#>  73719.81 279962.80 961517.23 761723.40

Los descriptivos apuntan a que los datos proporcionados en cada experimento tienen una escala muy distinta entre ellos. Estas diferencias se observan dentro de cada grupo experimental y entre ellos. Esto dificultará su comparación. De modo que habrá que normalizarlos o estandarizarlos para que puerdan ser comparados.

3.2. Representamos gráficamente los datos:

Previamente preparamos el input que precisa la función de representación del gráfico:

library(reshape)
mdata <- melt(datos, id=c("experiment","group"))
colnames(mdata) <- c("experiment", "group", "time", "growth")
dim(mdata)
#> [1] 68  4
attach(mdata)

Vemos cómo hemos organizado los datos:

head(mdata)
#>   experiment group time growth
#> 1       NE61    WT  t24 110000
#> 2       NE61    KO  t24 115000
#> 3       NE62    WT  t24 155000
#> 4       NE62    KO  t24 125000
#> 5       NE63    WT  t24 185000
#> 6       NE63    KO  t24 190000

El gráfico del crecimiento celular muestra patrones similares en ambos grupos, aunque las medianas (líneas gruesas dentro de cada caja) son mayores en todos los tiempos:

par(mfrow = c(1, 2))
plot(growth ~ time, col = "red",
     varwidth = TRUE, subset = group == "WT", main = "WT", ylim =c(0,3000000))
plot(growth ~ time, col = "green",
     varwidth = TRUE, subset = group == "KO", main = "KO", ylim =c(0,3000000))

plot of chunk unnamed-chunk-8

par(mfrow = c(1, 1))

3.3. Transformación logarítmica de los datos.

Una transformación logarítmica, permite visualizar mejor los datos cuando existen grandes diferencias de escala. Tras la aplicación de log2 sobre los datos, los representamos y se sigue observando el mismo patrón en ambos grupos pero con un incremento del crecimiento celular en KO:

mdata[, "log2_growth"]  <- log2(mdata$growth)
attach(mdata)
#> The following objects are masked from mdata (position 3):
#> 
#>     experiment, group, growth, time
par(mfrow = c(1, 2))
plot(log2_growth ~ time, col = "red",
     varwidth = TRUE, subset = group == "WT", main = "WT", ylim =c(12,22))
plot(log2_growth ~ time, col = "green",
     varwidth = TRUE, subset = group == "KO", main = "KO", ylim =c(12,22))

plot of chunk unnamed-chunk-9

par(mfrow = c(1, 1))

3.4. Normalización o estandarización de los datos.

Existen diferentes formas de normalizar los datos para que puedan ser comparados entre ellos. Una de ellas sería el cambio de escala en cada experimento, considerando como referencia el nivel de crecimiento en el tiempo 24 (referencia basal).

A continuación dividimos todos los datos de cada experimento por el valor de crecimiento en el tiempo 24:

ndatos  <- datos[,c("t24", "t48", "t72", "t96")] /  datos[,"t24"] 
ndatos <- cbind(datos[,c("experiment", "group")],ndatos)

Estos son los datos normalizados considerando como referencia basal el tiempo 24 en cada experimento:

ndatos
#>    experiment group t24       t48       t72       t96
#> 1        NE61    WT   1  2.272727  3.227273  7.909091
#> 2        NE61    KO   1  1.608696  3.884052  9.434783
#> 3        NE62    WT   1  1.838710  1.677419  3.741935
#> 4        NE62    KO   1  1.960000  2.320000  4.800000
#> 5        NE63    WT   1  8.324324 13.216216 10.918919
#> 6        NE63    KO   1  3.210526 10.473684 12.157895
#> 7        NE63    KO   1  7.318182 24.363636 13.636364
#> 8        NE64    WT   1  0.637931  1.448276  1.655172
#> 9        NE64    KO   1  1.382353  2.411765  2.784312
#> 10       NE78    WT   1  1.729730  2.243243  3.891892
#> 11       NE78    KO   1  1.218182  1.454545  2.745455
#> 12       NE79    WT   1  6.428571  8.071429 11.714286
#> 13       NE79    KO   1 10.233333 16.466667 38.800000
#> 14       NE80    WT   1  3.270270  8.364865 10.810811
#> 15       NE80    WT   1  3.756757 14.378378 20.540541
#> 16       NE80    KO   1  3.722222 20.111111 26.222222
#> 17       NE80    KO   1  3.318182 18.545455 25.454545

Describimos los datos normalizados:

mdata <- melt(ndatos, id=c("experiment","group"))
colnames(mdata) <- c("experiment", "group", "time", "growth")
attach(mdata)
#> The following objects are masked from mdata (position 3):
#> 
#>     experiment, group, growth, time
#> The following objects are masked from mdata (position 4):
#> 
#>     experiment, group, growth, time
summary(mdata)
#>  experiment group    time        growth       
#>  NE61: 8    KO:36   t24:17   Min.   : 0.6379  
#>  NE62: 8    WT:32   t48:17   1st Qu.: 1.0000  
#>  NE63:12            t72:17   Median : 2.9974  
#>  NE64: 8            t96:17   Mean   : 6.4574  
#>  NE78: 8                     3rd Qu.: 9.6344  
#>  NE79: 8                     Max.   :38.8000  
#>  NE80:16

Representamos gráficamente los datos normalizados y las diferencias se aprecian con mayor claridad:

par(mfrow = c(1, 2))
plot(growth ~ time, col = "red",
     varwidth = TRUE, subset = group == "WT", main = "WT", ylim =c(0,40))
plot(growth ~ time, col = "green",
     varwidth = TRUE, subset = group == "KO", main = "KO", ylim =c(0,40))

plot of chunk unnamed-chunk-13

par(mfrow = c(1, 1))



4. Detección de diferencias significativas: ANOVA de Friedman

La muestra presenta diferencias entre WT y KO a lo largo del tiempo. ¿Son diferencias estadísticamente significativas?. Aplicaremos el test de ANOVA de Friedman. Es una prueba no paramétrica que permite abordar este tipo de estudios cuando no podemos garantizar la normalidad de los datos en los distintos grupos evaluados:

resp <- aggregate(growth,
                  by = list(g = group,
                            t = time),
                  FUN = mean)
resp
#>    g   t         x
#> 1 KO t24  1.000000
#> 2 WT t24  1.000000
#> 3 KO t48  3.774631
#> 4 WT t48  3.532378
#> 5 KO t72 11.114546
#> 6 WT t72  6.578387
#> 7 KO t96 15.115064
#> 8 WT t96  8.897831
friedman.test(x ~ g | t, data = resp) 
#> 
#>  Friedman rank sum test
#> 
#> data:  x and g and t
#> Friedman chi-squared = 3, df = 1, p-value = 0.08326




5. Bibliografía y enlaces de interés: