¡Descarga Lección 24 R y más Apuntes en PDF de Matemáticas solo en Docsity!
Lección 24
ANOVA básico
En esta lección explicamos cómo efectuar con R los análisis de la varianza básicos que se han estudiado en el curso. Aunque el tema central es el estudio de los aspectos técnicos del ANOVA con R y los tests posteriores de comparación de pares de medias, dedicamos también una pequeña sección a introducir algunas funciones que sirven para comprobar dos de las condiciones del ANOVA —la normalidad y la igualdad de varianzas— sin pro- fundizar en las mismas. En la próxima lección estudiaremos con más detalle algunos tests de normalidad.
24.1. Los modelos del ANOVA en R
Los modelos a los que se aplica un ANOVA u otras muchas funciones, como por ejemplo la función lm para calcular la recta de regresión lineal que veíamos en la Lección 3, se especifican en R mediante fórmulas. El operador básico para construir una fórmula es la tilde ~. Las fórmulas suelen tener la forma Y~modelo, donde la Y es un vector y el modelo es una combinación de vectores que representa el modelo con el que queremos explicar el vector Y (para ser precisos, al que queremos ajustar el vector Y). Por ejemplo, cuando en la Lección 3 queríamos calcular la recta de regresión por mínimos cuadrados de un vector Y respecto de un vector X, usábamos lm(Y~X). Esto significa que aplicábamos la función lm a la fórmula Y~X que indica que queremos explicar Y como una función lineal de X.
En el ANOVA básico que hemos estudiado en el curso, usamos cuatro tipos de fórmulas. Sean X una variable numérica y F1 y F2 dos factores de una cierta tabla de datos.
La fórmula X~F1 se usa para indicar el ANOVA de un factor, F1, de la variable X.
La fórmula X~F1+F2 se usa para indicar el ANOVA de dos factores, F1 y F2, de la variable X, sin tener en cuenta la interacción entre los factores. Es decir, suponiendo que sus efectos se acumulan, sin que haya interacción entre los mismos.
La fórmula X~F1*F2 se usa para indicar el ANOVA de dos factores, F1 y F2, de la variable X, teniendo en cuenta además la interacción entre los dos factores.
La fórmula X~F1:F2 se usa para indicar el ANOVA de un factor que tiene como niveles los pares de niveles de F1 y F2.
La función básica de R para realizar un ANOVA es aov. Su sintaxis genérica es
aov(formula,data=...)
con los argumentos siguientes:
formula: Una fórmula que especifique un modelo de ANOVA.
data: Opcional, sirve para especificar, si es necesario, el data frame al que pertenecen las variables utilizadas en la fórmula.
Así, por ejemplo, si tenemos un data frame llamado DF, con una variable numérica X y un factor Fact, para realizar el ANOVA de la variable X respecto del factor Fact con la función aov podríamos entrar
aov(X˜Fact,data=DF) o aov(DF$X˜DF$Fact).
Otra posibilidad, que por ahora no usaremos pero que será útil más adelante, es aplicar la función anova (no confundirla con aov) al resultado de una función que calcule algún tipo de ajuste de una variable a un modelo, por ejemplo a una lm. La sintaxis sería entonces
anova(lm(formula,data=...))
24.2. ANOVA de un factor
Para ilustrar el ANOVA de un factor con R utilizaremos un experimento en el que se quiso determinar si cuatro dietas determinadas tenían alguna influencia en el tiempo de coagulación de la sangre en mamíferos.^1 Para ello se escogieron 24 animales, se repartieron de manera aleatoria en 4 grupos de 6 ejemplares cada uno, y a cada grupo se le asignó de manera aleatoria una de las 4 dietas bajo estudio. Al cabo de un cierto tiempo se midió el tiempo de coagulación de la sangre en estos animales. Los resultados se muestran en la Tabla 24.1. Observad que estamos ante un diseño experimental de un solo factor (las dietas) y de efectos fijos, ya que los niveles del factor (las dietas usadas) han sido escogidos por el experimentador.
Para determinar si los tiempos medios de coagulación son los mismos para las cuatro dietas o no, vamos a realizar un ANOVA sobre estos datos. Para ello, vamos a recoger estos datos en un data frame. (^1) Véase Statistics for Experimenters (2a edición), de G. P. Box, W. G. Hunter y J. S. Hunte (Wiley, 2005), p. 133.
Obtenemos el gráfico de la Figura 24.1, donde vemos que las medias para las dietas A, B y C tienen el aspecto de ser diferentes dos a dos, y en cambio es posible que no podamos rechazar que las de las dietas A y D sean iguales. Por tanto el resultado que esperamos del ANOVA es que nos permita rechazar la hipótesis nula de que los cuatro tiempos medios de coagulación son iguales.
A B C D
60
65
70
Figura 24.1. Diagrama de cajas de los tiempos de coagula- ción según las diferentes dietas.
Como hemos comentado, para realizar el ANOVA deseado, la instrucción básica es la siguiente:
aov(coag~diet,data=coagulacion) Call: aov(formula = coag ~ diet, data = coagulacion)
Terms: diet Residuals Sum of Squares 228 112 Deg. of Freedom 3 20
Residual standard error: 2. Estimated effects may be unbalanced
El resultado no es la tabla del ANOVA. Para obtenerla, hay que aplicar summary al resultado de aov:
summary(aov(coag~diet,data=coagulacion))
Df Sum Sq Mean Sq F value Pr(>F) diet 3 228 76.0 13.571 4.658e-05 *** Residuals 20 112 5.
Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1
El resultado de esta función summary es la tabla ANOVA usual:
En la primera columna, dos etiquetas: el nombre del factor, en este caso diet, y Residuals, que representa los errores o residuos del ANOVA.
La segunda columna, etiquetada Df, nos da los grados de libertad correspondientes al factor (su número de niveles menos 1) y a los residuos (el número de individuos en la tabla, menos el número de niveles del factor).
La tercera columna, etiquetada con Sum Sq, nos muestra las sumas de los cuadrados del factor, SST r, y de los residuos, SSE.
La cuarta columna, etiquetada con Mean Sq, contiene las medias de los cuadrados del factor, M ST r, y de los residuos, M SE.
La quinta columna, etiquetada con F value, nos da el valor del estadístico de con- traste.
En la sexta columna, etiquetada con Pr(>F), aparece el p-valor del contraste.
La séptima columna, sin etiqueta, indica el nivel de significación del p-valor según el código usual, explicado en la última línea. A mayor número de asteriscos, más significativo es el p-valor y por lo tanto hay más evidencia para rechazar la hipótesis nula.
En nuestro caso, hemos obtenido un p-valor del orden de 5 · 10 −^5. Esto nos permite rechazar la hipótesis nula y concluir que no todos los tiempos medios de coagulación para las diferentes dietas consideradas son iguales, como esperábamos. Si ahora queremos de- terminar qué pares de dietas dan tiempos medios de coagulación diferentes, tendremos que realizar algún test de comparaciones de pares de medias: véase la Sección 24.6.
Como hemos comentado, una manera alternativa de realizar un ANOVA es mediante anova(lm(...)). Con esta construcción obtenemos directamente la tabla, sin necesidad de aplicar summary.
anova(lm(coag~diet,data=coagulacion)) Analysis of Variance Table
Response: coag Df Sum Sq Mean Sq F value Pr(>F) diet 3 228 76.0 13.571 4.658e-05 ***
$ Mean Sq: num 76 5. $ F value: num 13.6 NA $ Pr(>F) : num 4.66e-05 NA
- attr(*, "heading")= chr "Analysis of Variance Table\n" "Response: coag"
tabla2$"F value" [1] 13.57143 NA tabla2$"Pr(>F)"[1] #el p-valor [1] 4.658471e-
Veamos otro ejemplo de ANOVA de un factor.
Ejemplo 24.1. En un experimento, se estudió el efecto de seis dietas sobre el crecimiento de crías de conejos domésticos.^2 Los datos obtenidos están recogidos en la tabla de datos rabbit del paquete faraway.
install.packages("faraway",dep=TRUE) #si aún no está instalado library(faraway) str(rabbit) ’data.frame’: 30 obs. of 3 variables: $ treat: Factor w/ 6 levels "a","b","c","d",..: 6 2 3 3 1 2 3 6 4 1 ... $ gain : num 42.2 32.6 35.2 40.9 40.1 38.1 34.6 34.3 37.5 44. ... $ block: Factor w/ 10 levels "b1","b10","b2",..: 1 1 1 3 3 3 4 4 4 5 ... help(rabbit)
La variable treat indica la dieta, sus valores son a, b,... , f , y la variable gain indica el aumento de peso; la variable block, que indica la camada, es irrelevante en este análisis concreto.
Si dibujamos el diagrama de cajas de los crecimientos para cada dieta,
boxplot(gain~treat,data=rabbit)
obtenemos el gráfico de la Figura 24.2, donde no observamos grandes diferencias en los aumentos medios de peso.
Para determinar si hay diferencia en los aumentos medios de peso bajo las seis dietas, realizaremos un ANOVA de la variable gain ajustándola al factor treat.
summary(aov(gain~treat,data=rabbit)) Df Sum Sq Mean Sq F value Pr(>F) treat 5 293.4 58.68 1.886 0. Residuals 24 746.5 31.
El p-valor es 0.134, lo que indica que, efectivamente, no hay evidencia que permita afirmar que los aumentos medios no son todos iguales. Si quisiéramos obtener sólo el p-valor, podríamos haber hecho simplemente
(^2) Véase Experimental Design and Analysis, de M. Lentner y T. Bishop (Valley Book Co. 1986), p. 428.
a b c d e f
20
25
30
35
40
45
50
Figura 24.2. Diagrama de cajas de los crecimientos según las diferentes dietas.
summary(aov(gain~treat,data=rabbit))[[1]]$"Pr(>F)"[1] [1] 0.
Al usar las instrucciones aov o anova(lm(...)) para realizar un ANOVA, se ha de emplear un factor (o varios, en las próximas secciones) para separar la variable numérica en subpoblaciones. Veamos un ejemplo de lo que pasa si nos descuidamos en este punto.
Ejemplo 24.2. En un experimento, se estudió el efecto de la vitamina C en el crecimiento de los dientes.^3 Se trató, en diferentes momentos, a cada ejemplar de un grupo de 10 cobayas con todas las combinaciones posibles de dosis de vitamina C (0.5, 1 o 2 mg) y métodos de suministro (mediante zumo de naranja o como ácido ascórbico) durante 6 semanas, y se cuantificó el crecimiento de sus dientes durante dicho período (más en concreto, se midió la longitud media de sus odontoblastos al final del mismo). El resultado es una tabla de 60 datos, que aparecen recogidos en el fichero ToothGrowth del paquete UsingR.
install.packages("UsingR",dep=TRUE) #si aún no está instalado library(UsingR) str(ToothGrowth) ’data.frame’: 60 obs. of 3 variables: $ len : num 4.2 11.5 7.3 5.8 6.4 10 11.2 11.2 5.2 7 ... $ supp: Factor w/ 2 levels "OJ","VC": 2 2 2 2 2 2 2 2 2 2 ... $ dose: num 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0. ... help(ToothGrowth)
(^3) Véase The Statistics of Bioassay de C. I. Bliss (Academic Press, 1952), p. 499–501.
había diferencias en su efectividad. Los cuatro procesos usaban una técnica de cultivo su- mergido y empleaban agua de macerado de maíz como fuente de nitrógeno orgánico. Como la composición de este líquido puede afectar la producción final de penicilina, para evaluar los cuatro procesos se prepararon 5 mezclas diferentes de agua de macerado de maíz, de cada mezcla se tomaron 4 muestras y se asignaron de manera aleatoria a los cuatro procesos de producción. Los resultados obtenidos son los de la Tabla 24.2, y aparecen recogidos en la tabla de datos penicillin del paquete faraway.
Proceso de prod. Mezcla A B C D 1 89 88 97 94 2 84 77 92 79 3 81 87 87 85 4 87 92 89 84 5 79 81 80 88
Tabla 24.2. Tabla de la producción de penicilina según el proceso de producción y la mezcla.
Observad que estamos ante un diseño experimental de bloques completos aleatorios. Los bloques son las mezclas de agua de macerado de maíz y los tratamientos (los procesos de producción) se han asignado de manera aleatoria a las unidades experimentales (las muestras) de cada bloque, de tal manera que cada bloque contiene exactamente una unidad experimental tratada con cada tratamiento. Este es un ejemplo paradigmático de uso de bloques: para evitar la influencia de la variable “extraña” dada por la composición del agua de macerado de maíz, que puede influir en la producción, se escogen algunas composiciones y se prueban todos los procesos de manera independiente sobre cada composición.
Demos un vistazo a esta tabla de datos.
str(penicillin) ’data.frame’: 20 obs. of 3 variables: $ treat: Factor w/ 4 levels "A","B","C","D": 1 2 3 4 1 2 3 4 1 2 ... $ blend: Factor w/ 5 levels "Blend1","Blend2",..: 1 1 1 1 2 2 2 2 3 3 ... $ yield: num 89 88 97 94 84 77 92 79 81 87 ... head(penicillin) treat blend yield 1 A Blend1 89 2 B Blend1 88 3 C Blend1 97 4 D Blend1 94 5 A Blend2 84 6 B Blend2 77
La variable treat contiene el proceso de producción, con valores A, B, C y D; la variable
blend contiene la mezcla usada (los bloques), con valores Blend1, Blend2,... , Blend5; y la variable numérica yield contiene un valor que cuantifica la producción de penicilina.
Ahora es un buen momento para mencionar que a partir de un data frame con dos factores podemos obtener su tabla de doble entrada, con filas y columnas representando los niveles de ambos factores (la Tabla 24.2, en este ejemplo), por medio de la instrucción
ftable(xtabs(formula,data=...))
Por ejemplo
ftable(xtabs(yield~blend+treat,data=penicillin)) treat A B C D blend Blend1 89 88 97 94 Blend2 84 77 92 79 Blend3 81 87 87 85 Blend4 87 92 89 84 Blend5 79 81 80 88
produce la Tabla 24.2. Observad el orden en el que hemos entrado los factores en la fórmula: el primero, en nuestro caso blend, indica las filas y el segundo, treat, las columnas.
Veamos como son los diagramas de cajas de la producción de penicilina separada por procesos de producción y por mezclas.
boxplot(yield~treat,data=penicillin) boxplot(yield~blend,data=penicillin)
Obtenemos los gráficos de la Figura 24.3. Vemos que no hay mucha diferencia entre las producciones para los diferentes procesos (sin tener en cuenta los bloques), y que sí que hay algunas diferencias en las producciones según la composición del agua de macerado de maíz.
Si realizamos un ANOVA de un factor para cada uno de estos dos factores por separado, obtenemos resultados consistentes con esta observación visual:
summary(aov(yield~treat,data=penicillin)) Df Sum Sq Mean Sq F value Pr(>F) treat 3 70 23.33 0.762 0. Residuals 16 490 30. summary(aov(yield~blend,data=penicillin)) Df Sum Sq Mean Sq F value Pr(>F) blend 4 264 66.00 3.345 0.038 * Residuals 15 296 19.
Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1
el factor a analizar, y los procesos de producción como los bloques. En efecto:
summary(aov(yield~blend+treat,data=penicillin)) Df Sum Sq Mean Sq F value Pr(>F) blend 4 264 66.00 3.504 0.0407 * treat 3 70 23.33 1.239 0. Residuals 12 226 18.
Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1
Podemos extraer los resultados de una tabla de un ANOVA de bloques generada por summary(aov(...)) añadiendo los mismos sufijos que en el caso de un factor:
tabla=summary(aov(yield~treat+blend,data=penicillin)) tabla Df Sum Sq Mean Sq F value Pr(>F) treat 3 70 23.33 1.239 0. blend 4 264 66.00 3.504 0.0407 * Residuals 12 226 18.
Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1
tabla[[1]]$"Mean Sq" [1] 23.33333 66.00000 18. tabla[[1]]$"Pr(>F)" [1] 0.33865812 0.04074617 NA tabla[[1]]$"Pr(>F)"[1] [1] 0.
Veamos otro ejemplo de ANOVA de bloques completos aleatorios.
Ejemplo 24.3. Para estudiar si las diferentes actividades de evaluación llevadas a cabo en la asignatura de Matemáticas I tienen una dificultad similar (o mejor dicho, para mirar de confirmar nuestras sospechas de que no es así), escogimos una muestra aleatoria simple de 15 estudiantes de Biología o Bioquímica que se hubieran presentado al control 2 de dicha asignatura en el curso 2012-2013, y anotamos las notas obtenidas por estos estudiantes en los apartados de Tests, Talleres, Ejercicios de Casa, Control 1 y Control 2. Si obtenemos evidencia de que no todas las medias de las notas de las diferentes actividades fueron iguales, podremos concluir que no todas las actividades tienen la misma dificultad. Los datos obtenidos son los de la tabla 24.3.
Como podemos observar, se trata de un experimento de bloques completos aleatorios: hemos escogido de manera aleatoria unos bloques (los estudiantes) y sobre cada estudiante hemos calculado el valor de cada uno de los niveles del factor a estudiar (las notas en las diferentes actividades). Este diseño es el adecuado para este problema, puesto que hay una gran variabilidad en las notas obtenidas por estudiantes diferentes, desde matrículas a
Actividad Estudiante Tests Casa Talleres Control 1 Control 2 1 48 42 85 31 20 2 94 86 100 52 48 3 98 94 93 93 90 4 60 58 71 66 46 5 52 56 79 64 24 6 77 73 84 80 95 7 84 84 94 83 70 8 83 76 99 51 55 9 75 55 70 85 48 10 24 49 57 19 55 11 47 47 64 44 35 12 53 42 78 50 8 13 82 93 92 66 30 14 66 62 74 24 56 15 49 40 74 35 21
Tabla 24.3. Tabla de las notas obtenidas por 15 estudiantes en las diferentes actividades de evaluación.
suspensos. Al tomar bloques, es decir, al considerar todas las notas de un conjunto aleatorio fijo de estudiantes, eliminamos el efecto de esta variabilidad.
Vamos a construir un data frame con estos datos. Este data frame tendrá tres variables: notas, con las notas obtenidas por los estudiantes; acts, con los tipos de actividades de evaluación realizados; y bloques, con el indicador de cada estudiante/bloque.
notas=c(48,42,85,31,20,94,86,100,52,48,98,94,93,93,90,60,58,71,66,46, 52,56,79,64,24,78,73,84,80,95,84,84,94,83,70,83,76,99,51,55, 75,55,70,85,48,24,49,57,19,55,47,47,64,44,35,53,42,78,50,8, 82,93,92,66,30,66,62,74,24,56,49,40,74,35,21) acts=rep(c("Tests","Casa","Talleres","Control1","Control2"),15) bloques=as.factor(rep(1:15,each=5)) #los bloques han de ser un factor notas.bloques=data.frame(notas,acts,bloques) str(notas.bloques) ’data.frame’: 75 obs. of 3 variables: $ notas : num 48 42 85 31 20 94 86 100 52 48 ... $ acts : Factor w/ 5 levels "Casa","Control1",..: 5 1 4 2 3 5 1 4 2 3 ... $ bloques: Factor w/ 15 levels "1","2","3","4",..: 1 1 1 1 1 2 2 2 2 2 ... head(notas.bloques) notas acts bloques 1 48 Tests 1 2 42 Casa 1
Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1
El p-valor que nos interesa es el de la primera fila, etiquetada con el factor cuyos niveles queremos comparar, en este caso acts. Este p-valor es muy pequeño, 1. 31 · 10 −^7 , lo que indica que, en efecto, no todas las notas medias de las actividades de evaluación fueron iguales.
24.4. ANOVA de dos vías
Para ilustrar el ANOVA de dos vías, es decir, de dos factores con interacción, analizaremos los resultados de un experimento sobre el efecto de venenos y antídotos.^5 En este experi- mento se usaron tres venenos y cuatro antídotos, y cada combinación de veneno y antídoto se administró a cuatro ratas elegidas de manera aleatoria e independiente. A continuación, se anotó el tiempo de supervivencia de cada animal, en unidades de 10 horas. Los resul- tados se muestran en la Tabla 24.4. El objetivo era analizar el efecto de los venenos y de los antídotos. Se trata, pues, de un diseño experimental de dos factores, en el que cada individuo ha sido asignado al azar a cada nivel de los dos factores. Estamos interesados en contrastar la igualdad de medias bajo estos dos factores.
Antídoto Veneno A B C D I 0.31 0.82 0.43 0. 0.45 1.10 0.45 0. 0.46 0.88 0.63 0. 0.43 0.72 0.76 0. II 0.36 0.92 0.44 0. 0.29 0.61 0.35 1. 0.40 0.49 0.31 0. 0.23 1.24 0.40 0. III 0.22 0.30 0.23 0. 0.21 0.37 0.25 0. 0.18 0.38 0.24 0. 0.23 0.29 0.22 0.
Tabla 24.4. Tabla de tiempos de supervivencia según el veneno y el antídoto.
Los datos de este experimento están contenidos en la tabla rats del paquete faraway. Vamos a explorar estos datos.
(^5) Véase “An analysis of transformations”, de G. Box y D. Cox, J. Roy. Stat. Soc. Series B 26 (1964), pp. 211–252. Véase también Practical Regression and Anova using R, de J. Faraway, p. 182.
str(rats) ’data.frame’: 48 obs. of 3 variables: $ time : num 0.31 0.82 0.43 0.45 0.45 1.1 0.45 0.71 0.46 0. ... $ poison: Factor w/ 3 levels "I","II","III": 1 1 1 1 1 1 1 1 1 1 ... $ treat : Factor w/ 4 levels "A","B","C","D": 1 2 3 4 1 2 3 4 1 2 ... head(rats) time poison treat 1 0.31 I A 2 0.82 I B 3 0.43 I C 4 0.45 I D 5 0.45 I A 6 1.10 I B help(rats)
La variable treat contiene el antídoto, con valores A, B, C y D; la variable poison contiene el veneno, con valores I, II y III; y la variable numérica time contiene el tiempo de supervivencia de las ratas. Veamos cómo son los diagramas de cajas de estos tiempos, separados por venenos y por antídotos.
boxplot(time~poison,data=rats) boxplot(time~treat,data=rats)
Obtenemos los gráficos de la Figura 24.5. Parece que hay diferencias en los tiempos de supervivencia tanto para los diferentes venenos como para los diferentes antídotos. Podemos también dibujar un diagrama de cajas de los tiempos de supervivencia separándolos por combinaciones de veneno y antídoto. La instrucción para hacerlo es la siguiente (observad la fórmula del argumento) y el resultado es el de la Figura 24.6.
boxplot(time~poison:treat,data=rats)
El código para llevar a cabo el ANOVA de dos vías del tiempo de supervivencia de las ratas bajo los efectos combinados de veneno y antídoto es el siguiente:
summary(aov(time~poison*treat,data=rats)) Df Sum Sq Mean Sq F value Pr(>F) poison 2 1.0330 0.5165 23.222 3.33e-07 *** treat 3 0.9212 0.3071 13.806 3.78e-06 *** poison:treat 6 0.2501 0.0417 1.874 0. Residuals 36 0.8007 0.
Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1
En la última columna de esta tabla, el primer p-valor es el del contraste de la igualdad de tiempos medios de supervivencia según el veneno (poison): el p-valor 3. 3 · 10 −^7 indica
I.A II.A III.A I.B II.B III.B I.C II.C III.C I.D II.D III.D
Figura 24.6. Diagramas de cajas de tiempos de superviven- cia según diferentes combinaciones de veneno y antídoto.
la variable X para cada nivel de F 2 en el nivel de F 1 correspondiente. Si no hay ninguna interacción entre estos factores, las líneas resultantes serán paralelas. Cuanto más se alejen de ser paralelas, más evidencia de interacción habrá entre estos dos factores. Estos gráficos de interacción se dibujan en R con la instrucción interaction.plot aplicada, por este orden, a los dos factores y la variable numérica.
Así, por ejemplo, para obtener el gráfico de interacción de la variable antídoto respecto de la variable veneno en el ANOVA anterior, entraríamos
interaction.plot(rats$treat,rats$poison,rats$time)
y obtendríamos el gráfico de la izquierda en la Figura 24.7. En él observamos que los tres niveles del veneno producen líneas casi paralelas, aunque se observa una ligera interacción: la pendiente de la recta correspondiente al veneno II entre los valores medios de C y D es mucho mayor que la de las rectas correspondientes a los otros dos venenos.
Este diagrama no ha quedado muy bonito, ya que en él aparecen los nombres enteros de las variables, y en la caja de información los factores aparecen ordenados por altura, no por nombre. Una posibilidad para mejorarlo es cambiar las etiquetas (la etiqueta del factor en la caja de información se define con trace.label) y reordenar los factores en dicha caja (con fixed=TRUE). Ya que estamos, dibujaremos algo más gruesas las líneas.
interaction.plot(rats$treat,rats$poison,rats$time,xlab="Antídoto", ylab="Tiempo medio de supervivencia",trace.label="Veneno", lwd=c(2,2,2),fixed=TRUE)
Obtenemos el gráfico de la derecha en la Figura 24.7.
rats$treat
mean of rats$time
A B C D
rats$poison II I III
Antídoto
Tiempo medio de supervivencia
A B C D
Veneno I II III
Figura 24.7. Gráficos de interacción del veneno y el antído- to.
Veamos otro ejemplo de ANOVA de dos vías.
Ejemplo 24.4. En el Ejemplo 24.2 hablamos de un cierto experimento sobre el efecto de la vitamina C en el crecimiento de los dientes, y de la tabla de datos ToothGrowth del paquete UsingR que recoge los datos obtenidos en ese experimento. Ahora vamos a realizar un ANOVA de dos vías sobre esta tabla de datos para contrastar la influencia en dicho crecimiento de la dosis de vitamina C y del método de suministrarla, teniendo en cuenta que sus efectos se pueden combinar. Recordad de aquel ejemplo que, antes de llevar a cabo un ANOVA sobre esta tabla que involucre la variable dose que contiene la dosis, hay que convertir esta variable en un factor. Para ello, copiamos la tabla ToothGrowth en una nueva tabla ToothGrowth2, y en esta última modificamos dicha variable.
ToothGrowth2=ToothGrowth ToothGrowth2$dose=as.factor(ToothGrowth2$dose) str(ToothGrowth2) ’data.frame’: 60 obs. of 3 variables: $ len : num 4.2 11.5 7.3 5.8 6.4 10 11.2 11.2 5.2 7 ... $ supp: Factor w/ 2 levels "OJ","VC": 2 2 2 2 2 2 2 2 2 2 ... $ dose: Factor w/ 3 levels "0.5","1","2": 1 1 1 1 1 1 1 1 1 1 ...
Demos ahora una ojeada a los diagramas de cajas de esta tabla, separando las longitudes de los odontoblastos por dosis, por suministro y por ambos.
boxplot(len~dose,data=ToothGrowth2)
