ANOVA

DISEÑO DE UN FACTOR

Comparar tres poblaciones

Iris.setosa <- c(5.1, 4.9, 4.7, 4.6, 5, 5.4, 4.6, 5, 4.4, 4.9,5.4, 4.8, 4.8, 4.3, 5.8)
Iris.versicolor <- c(7, 6.4, 6.9, 5.5, 6.5, 5.7, 6.3, 4.9, 6.6,5.2, 5, 5.9, 6, 6.1, 5.6)
Iris.virginica <- c(6.3, 5.8, 7.1, 6.3, 6.5, 7.6, 4.9, 7.3, 6.7,7.2, 6.5, 6.4, 6.8, 5.7, 5.8)

Sin embargo con los datos anterioes no se puede trabajar por lo tanto se realiza la union de los tres grupos en un solo vector añadiendo el nombre de las variables(factores)

longitud <- c(Iris.setosa, Iris.versicolor, Iris.virginica)#unión de los 3 vectores anterioes
especie <- rep(1:3, each = 15)#crear los 3 variable las cuales se repitan 15 veces 
especie <- factor(especie, labels = c("Iris setosa", "Iris versicolor","Iris virginica"))#Agregar nombre de las 3 variables por cada vector
especie <- gl(3, 15, labels = c("Iris setosa", "Iris versicolor","Iris virginica"))#definir como factor para poder trabajarlos
split(longitud, especie)#Ver los vectores por variable

## $`Iris setosa`
##  [1] 5.1 4.9 4.7 4.6 5.0 5.4 4.6 5.0 4.4 4.9 5.4 4.8 4.8 4.3 5.8
## 
## $`Iris versicolor`
##  [1] 7.0 6.4 6.9 5.5 6.5 5.7 6.3 4.9 6.6 5.2 5.0 5.9 6.0 6.1 5.6
## 
## $`Iris virginica`
##  [1] 6.3 5.8 7.1 6.3 6.5 7.6 4.9 7.3 6.7 7.2 6.5 6.4 6.8 5.7 5.8

tapply(longitud, especie, summary)#Resumen

## $`Iris setosa`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   4.300   4.650   4.900   4.913   5.050   5.800 
## 
## $`Iris versicolor`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   4.900   5.550   6.000   5.973   6.450   7.000 
## 
## $`Iris virginica`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    4.90    6.05    6.50    6.46    6.95    7.60

plot(longitud ~ especie)#Gráfico de cajas(boxplot)

Asumiendo que la variable longitud sigue una distribución normal con varianza común para las tres poblaciones, el procedimiento es:

1. POR ANáLISIS DE VARIANZA(aov)

p.aov <- aov(longitud ~ especie)#comparación de la lonitud por especie(la especie es nuestro unico factor la cual tiene 3 niveles(Iris setosa,Iris versicolor y Iris virginica))
summary(p.aov)#resumen

##             Df Sum Sq Mean Sq F value  Pr(>F)    
## especie      2  18.76   9.382   25.71 5.1e-08 ***
## Residuals   42  15.32   0.365                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model.tables(p.aov)#ver el valor p por cada variable

## Tables of effects
## 
##  especie 
## especie
##     Iris setosa Iris versicolor  Iris virginica 
##         -0.8689          0.1911          0.6778

model.tables(p.aov, type = "mean")#ver la media general y por variable 

## Tables of means
## Grand mean
##          
## 5.782222 
## 
##  especie 
## especie
##     Iris setosa Iris versicolor  Iris virginica 
##           4.913           5.973           6.460

2. POR MODELO LINEAL (lm)

Otra posibilidad es definir el modelo lineal y obtener la tabla con la instrucci¶on anova.

g.lm <- lm(longitud ~ especie)#formula para comparación de la lonitud por especie
anova(g.lm)#mostrar anova

## Analysis of Variance Table
## 
## Response: longitud
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## especie    2 18.763  9.3816  25.715 5.105e-08 ***
## Residuals 42 15.323  0.3648                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(g.lm)#mostrar el resumen

## 
## Call:
## lm(formula = longitud ~ especie)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.56000 -0.31333 -0.01333  0.42667  1.14000 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              4.9133     0.1560  31.505  < 2e-16 ***
## especieIris versicolor   1.0600     0.2206   4.806 1.99e-05 ***
## especieIris virginica    1.5467     0.2206   7.013 1.39e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.604 on 42 degrees of freedom
## Multiple R-squared:  0.5505, Adjusted R-squared:  0.5291 
## F-statistic: 25.72 on 2 and 42 DF,  p-value: 5.105e-08

Por ambos procedimientos llegamos a los mismos resultados, sin embargo el modelo lineal contiene mayor información con el summary #Como el p-valor es muy pequeño se concluye que hay diferencias muy significativas entre las tres especie, es decir no son iguales entre sí Despues se puede obtenr el error cuadrático medio o estimación insesgada de la varianza del modelo

ECM <- deviance(p.aov)/p.aov$df.residual#ECM por el oav
summary(g.lm)$sigma^2#ECM por el metodo lineal

## [1] 0.3648254

EJEMPLO

Análisis de la varianza con un único factor tratamiento y cuatro niveles (P,A,B,AB).

P <- c(10, 0, 15, -20, 0, 15, -5, NA, NA, NA)#VECTOR A
A <- c(20, 25, 33, 25, 30, 18, 27, 0, 35, 20)#VECTOR B
B <- c(15, 10, 25, 30, 15, 35, 25, 22, 11, 25)#VECTOR C
AB <- c(10, 5, -5, 15, 20, 20, 0, 10, NA, NA)#VECTOR D
descenso <- c(P, A, B, AB)#Unión de vectores
tratam <- gl(4, 10, labels = c("placebo", "fármaco A", "fármaco B","asociación AB"))#agregar variables

#Suponiendo normalidad y homogeneidad de las varianzas, planteamos el test sobre la igualdad de medias.
#Un resumen numérico y gráfico se puede obtener con las instrucciones
mean(descenso, na.rm = TRUE)#media general, sin considerar NA

## [1] 15.31429

tapply(descenso, tratam, summary)#Resumen po cada fármaco

## $placebo
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
## -20.000  -2.500   0.000   2.143  12.500  15.000       3 
## 
## $`fármaco A`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00   20.00   25.00   23.30   29.25   35.00 
## 
## $`fármaco B`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    10.0    15.0    23.5    21.3    25.0    35.0 
## 
## $`asociación AB`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##  -5.000   3.750  10.000   9.375  16.250  20.000       2

stripchart(descenso ~ tratam, method = "stack")#gráficp

#El modelo lineal y la tabla del análisis de la varianza son
g.lm <- lm(descenso ~ tratam)
anova(g.lm)

## Analysis of Variance Table
## 
## Response: descenso
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## tratam     3 2492.6  830.87  8.5262 0.0002823 ***
## Residuals 31 3020.9   97.45                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#El p-valor es inferior (0.0002823) al nivel de significación propuesto (0.01) de modo que rechazamos la hipótesis nula de igualdad de medias y admitimos que hay diferencias entre los fármacos.

#Para ver si hay diferencias entre los fármacos A y B calcularemos el intervalo de confianza para la diferencia de medias:
mediaA <- mean(descenso[tratam == "f¶armaco A"])#media del fármaco A
mediaB <- mean(descenso[tratam == "f¶armaco B"])#media del fármaco B
dif <- mediaA - mediaB#diferencia de medias entre farmaco A Y B
ee.dif <- summary(g.lm)$sigma * sqrt(1/10 + 1/10)
c(dif - qt(0.995, 31) * ee.dif, dif + qt(0.995, 31) * ee.dif)

## [1] NaN NaN

#Como este intervalo contiene al cero, podemos pensar que las diferencias entre A y B no son significativas. Aunque puede comprobarse que ambos fármacos difieren significativamente del placebo, cuando se realiza más de una comparación necesitamos un método de comparaciones múltiples (método de la diferencia significativa honesta de Tukey con la función TukeyHSD).

DISEÑO DE DOS FACTORES SIN INTERACCIÓN Se trata de un dise~no de bloques aleatorizados (cada ¯nca es un bloque).

produc <- c(2.1, 2.2, 1.8, 2, 1.9, 2.2, 2.6, 2.7, 2.5, 2.8, 1.8,1.9, 1.6, 2, 1.9, 2.1, 2, 2.2, 2.4, 2.1)#Vector de producción
fert <- gl(4, 5)#vector de 4 fertilizantes las cuales cada una se repite 5 veces
finca <- factor(rep(1:5, 4))#crear una vector de 1 al 5 que se repita 4 veces, las cuales seran las filas para lo siguiente
xtabs(produc ~ finca + fert)#Crear una matriz 4 columnas por 5 filas

##      fert
## finca   1   2   3   4
##     1 2.1 2.2 1.8 2.1
##     2 2.2 2.6 1.9 2.0
##     3 1.8 2.7 1.6 2.2
##     4 2.0 2.5 2.0 2.4
##     5 1.9 2.8 1.9 2.1

Ahora podemos generar un resumen de los datos y los gráficos

tapply(produc, fert, summary)#Resumen por variable(columna, fertilizante)

## $`1`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     1.8     1.9     2.0     2.0     2.1     2.2 
## 
## $`2`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2.20    2.50    2.60    2.56    2.70    2.80 
## 
## $`3`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.60    1.80    1.90    1.84    1.90    2.00 
## 
## $`4`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2.00    2.10    2.10    2.16    2.20    2.40

tapply(produc, finca, summary)#Resumen por variable(fila, finca)

## $`1`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.800   2.025   2.100   2.050   2.125   2.200 
## 
## $`2`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.900   1.975   2.100   2.175   2.300   2.600 
## 
## $`3`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.600   1.750   2.000   2.075   2.325   2.700 
## 
## $`4`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   2.000   2.000   2.200   2.225   2.425   2.500 
## 
## $`5`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.900   1.900   2.000   2.175   2.275   2.800

stripchart(produc ~ fert, method = "stack")#Gráfico de fertilizantes

stripchart(produc ~ finca, method = "stack")#Gráfico de fincas

interaction.plot(fert, finca, produc, legend = F)#interacción por fincas

interaction.plot(finca, fert, produc, legend = F)#interacción por fertilizantes

A la vista de los gráficos, concluimos que no hay datos atípicos, asimetrías o heterocedasticidad. Tampoco parece haber interacciones.

Modelo lineal y tabla de varianza son:

g.lm <- lm(produc ~ finca + fert)#los 2 factores son: finca y fertilizante
anova(g.lm)#Tabla de Varianza, la cual muestra que con un grado de signifcancia de 0.05 unicamente hay diferencia significativa entre fertilizante(0.00031) la cual tiene el valor p<0.05, mientras que entre finca no hay diferencias ya que el valor p>0.05 al tener un valor de 0.6395

## Analysis of Variance Table
## 
## Response: produc
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## finca      4  0.088 0.02200  0.6471 0.6395716    
## fert       3  1.432 0.47733 14.0392 0.0003137 ***
## Residuals 12  0.408 0.03400                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Los efectos y las medias son:

p.aov <- aov(produc ~ finca + fert)
efectos <- model.tables(p.aov)
medias <- model.tables(p.aov, type = "means")

En este caso el modelo es balanceado, de forma que el diseño es ortogonal y el orden de los factores en la instrucción anova no es importante. En este sentido hay que señalar que la tabla ANOVA de R corresponde a un contraste secuencial de modelos: y ~ 1 y ~ finca y ~ finca + fert

DISEÑO DE DOS FACTORES CON INTERACCIÓN

huevos <- c(93, 94, 93, 90, 93, 86, 95.5, 83.5, 92, 92.5, 82,82.5, 92, 91, 90, 95, 84, 78, 83.3, 87.6, 81.9, 80.1, 79.6,49.4, 84, 84.4, 77, 67, 69.1, 88.4, 85.3, 89.4, 85.4, 87.4,52, 77)#Vector numerico
genotipo <- rep(rep(1:3, each = 6), 2)#primer factor en vector numerico llamado genotipo , que tiene tres grupos repetidos 6 veces y estos 3 grupos repetidos 2 veces
siembra <- rep(1:2, each = 18)#Segundo factor en vector numerico llamado siembra que tiene 2 grupos repetido 8 veces
genotipo <- factor(genotipo, labels = c("++", "+-", "--"))#Tercer factor la cual remplaza los caracteres genotipicas al primer vector numerico del mismo nombre 
siembra <- factor(siembra, labels = c("100", "800"))#Segundo factor la cual es remplazada con los valores de siembra a emplear

El número de huevos eclosionados por casilla sigue la distribución binomial con n = 100 o n = 800. Para normalizar la muestra se aplica la transformación

y <- asin(sqrt(huevos/100))
y <- y * 180/pi
#de donde resulta la tabla:
split(round(y, 2), genotipo)

## $`++`
##  [1] 74.66 75.82 74.66 71.57 74.66 68.03 65.88 69.38 64.82 63.51 63.15
## [12] 44.66
## 
## $`+-`
##  [1] 77.75 66.03 73.57 74.11 64.90 65.27 66.42 66.74 61.34 54.94 56.23
## [12] 70.09
## 
## $`--`
##  [1] 73.57 72.54 71.57 77.08 66.42 62.03 67.46 71.00 67.54 69.21 46.15
## [12] 61.34

Aunque no es absolutamente necesario, vamos a poner los datos en forma de data.frame o base de datos de R.

problema <- data.frame(y, siembra, genotipo)#union de vectores y convercion a data.frame
rm(y, siembra, genotipo)#Remover los 3 objetos despues de convertirlos en data.fram
attach(problema)

gráficos

boxplot(y ~ siembra)#Diagrama de cajas comparando el factor siembra (cuando las cajas no se sonbreponen hay diferencias significativas)

boxplot(y ~ genotipo)#Diagrama de cajas comparando el factor genotipo (cuando las cajas se sonbreponen no hay diferencias significativas)

plot.design(problema, fun = "mean")#diagrama que indica la media de los factores

plot.design(problema, fun = "median")#diagrama que indica la mediana de los factores

interaction.plot(genotipo, siembra, y)#Gráfico que muestra la interacción genotip x siembra (cuando las lineas no son paralelas hay diferencias significativas)

interaction.plot(siembra, genotipo, y)#Gráfico que muestra la interacción siembra x genotipo(cuando las lineas son casi paralelas hay diferencias significativas)

Tabla del análisis de la varianza

p.aov <- aov(y ~ siembra * genotipo, data = problema)
summary(p.aov)

##                  Df Sum Sq Mean Sq F value   Pr(>F)    
## siembra           1  662.1   662.1  14.833 0.000574 ***
## genotipo          2    7.7     3.8   0.086 0.917952    
## siembra:genotipo  2   35.4    17.7   0.396 0.676456    
## Residuals        30 1339.1    44.6                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

No es significativa la diferencia entre los genotipos ni la interacción, pero sí existen diferencias significativas sembrando 100 u 800 huevos, siendo el porcentaje de eclosiones mayor en el primer caso (al haber menos huevos, las larvas disponen de más alimento). Esto se alcanza a ver en los boxplot e interaccion.plot

medias <- model.tables(p.aov, type = "means")#presentación de las medias
medias$tables$siembra

## siembra
##      100      800 
## 71.34575 62.76873

detach(problema)

DISEÑO DE CUADROS LATINOS

produc <- c(12, 17, 24, 12, 18, 22, 14, 15, 15, 13, 20, 31, 20,14, 12, 18)#Datos
fila <- factor(rep(1:4, 4))#Crear 4 filas que repitan 4 veces
columna <- factor(rep(1:4, each = 4))#Crear 4 columnas que repitan 4 veces
variedad <- c("A", "C", "D", "B", "B", "D", "C", "A", "C", "A","B", "D", "D", "B", "A", "C")#Datos de las variables
problema <- data.frame(fila, columna, variedad, produc)#Unión de los  4 factores
rm(fila, columna, variedad, produc)#remover los primeros factores que presentan vectores y quedarnos con la union de estas  que es con la que se va a trabajar
attach(problema)

matrix(problema$variedad, 4, 4)#agregar una matriz de 4x4 con las los variadades a trabajar para ver el diseño de cuadros latinos

##      [,1] [,2] [,3] [,4]
## [1,] "A"  "B"  "C"  "D" 
## [2,] "C"  "D"  "A"  "B" 
## [3,] "D"  "C"  "B"  "A" 
## [4,] "B"  "A"  "D"  "C"

La tabla del análisis de la varianza para este diseño es

p.aov <- aov(produc ~ fila + columna + variedad, data = problema)
summary(p.aov)

##             Df Sum Sq Mean Sq F value Pr(>F)  
## fila         3  18.69    6.23   0.527 0.6796  
## columna      3  35.19   11.73   0.993 0.4574  
## variedad     3 280.69   93.56   7.921 0.0165 *
## Residuals    6  70.88   11.81                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

graficos

boxplot(produc~fila)

boxplot(produc~columna)

boxplot(produc~variedad)

No hay diferencias signifícativas entre filas ni entre columnas. En cambio sí hay diferencias entre variedades.

DISEÑO MULTIFACTORIAL