Diseño de Experimentos. Guía docente.

Unidades 4 y 5

Introducción
Esta guía cubre las unidades 4 y 5. Estudiaremos principalmente como verificar el cumplimiento de las
suposiciones básicas del ANOVA y el cálculo del tamaño muestral.

Unidad 4. Comprobación de las suposiciones del modelo

En resumen, mediante el ANOVA analizamos los datos experimentales comparando las respuestas medias
en los distintos tratamientos y en caso de detectar diferencias significativas procedemos a un análisis más
detallado basado en comparaciones por parejas u otros contrastes. Todos estos procedimientos estadísticos se
basan en la suposición que los datos se ajustan al modelo

yij = µ + αi + eij

donde µ, αi son parámetros (cantidades fijas pero desconocidas) y eij son variables aleatorias independientes
y distribuidas según una normal de media 0 y varianza constante σ 2 . Por tanto, será necesario chequear si
estas suposiciones se satisfacen razonablemente.
Las tres suposiciones (o premisas) básicas que debemos validar son que los errores eij son:
1. independientes,
2. distribuidos según una normal,
3. la varianza es constante (no varía a través de los tratamientos)
La suposición de independencia es la más importante y difícil de corregir en caso que no se cumpla. Es
el diseño del experimento, un muestreo representativo, una adecuada aleatorización en la asignación de
los tratamientos y, en fin, cada paso del experimento, el que debe asegurarnos que las observaciones son
independientes. En esta unidad vamos a centrarnos en la validación de las otras dos premisas: normalidad y
varianza constante.
Veremos que la suposición de normalidad, aunque importante, puede relajarse en caso de muestras grandes.
La condición de varianza constante (también llamada condición de homocedasticidad) puede afectar sus-
tancialmente nuestras inferencias caso de no cumplirse. No obstante, la falta de homocedasticidad puede
corregirse en muchas situaciones. En las diversas secciones de este tema se consideran distintos métodos de
diagnóstico y posibles correcciones cuando las premisas del modelo no se cumplen.
Hasta cierto punto, podría parecer que “forzamos la máquina” porque vamos a estudiar una serie de
transformaciones con el fin que nuestros modelos lineales, con errores distribuidos normalmente, sean
aplicables a muchos tipos de datos. Cabe decir, que hay otras metodologías que podríamos usar en su lugar,
por ejemplo: los modelos lineales generalizados, los métodos robustos, los métodos de permutaciones o los
métodos no paramétricos basados en rangos que brevemente introduciremos al final del bloque. De hecho,
para ciertos tipos de datos, algunos de estos métodos alternativos pueden ser considerablemente más eficientes
(por ejemplo, producir intervalos de confianza más cortos con el mismo nivel de confianza) que el modelo
lineal basado en la distribución normal. Sin embargo, cada uno de estos métodos alternativos que acabamos
de mencionar podría ser motivo de un curso específico.

1
Evaluando el cumplimiento de los supuestos de modelización
Estrictamente las suposiciones de independencia, distribución normal de los errores con varianza constante
quizá no se cumplan por completo en datos reales. No obstante, los modelos lineales producen buenas
inferencias siempre que la desviación de los supuestos no sea demasiado grande. Por este motivo podemos
utilizar métodos gráficos y descriptivos para evaluar el grado de cumplimiento de las suposiciones. Como
veremos en esta unidad, también disponemos de contrastes estadísticos a tal fin.
Evaluamos el cumplimiento de las suposiciones sobre los errores aleatorios basándonos en los residuos. Los
residuos rij vienen dados por la diferencia entre el valor observado yij y el valor ajustado según el modelo
para dicha observación. En el caso del ANOVA de un factor, el valor ajustado para la observación yij es ȳi. .
Por tanto, para cada dato obtenemos un residuo

rij = yij − ȳi.

Los residuos rij son una muestra de la variable aleatoria eij . Según las suposiciones del modelo se asume
que eij ∼ N (0, σ 2 ), por tanto es a partir del estudio de los residuos que podremos validar los supuestos del
modelo.

Ejemplo
Se basa en el enunciado del problema 3.3 en la página 62 del libro de Oehlert que ya hemos utilizado en la
unidad 2.
Se trata de una diseño de un factor (factor Number of passes) con 5 niveles o tratamientos (0, 25, 75, 200,
500). La variable respuesta es la altura promedio (cm) de la vegetación a lo largo del trayecto un año después
del inicio del experimento. Un total de 20 trayectos (0.5 m ancho x 1.5 m largo) se asignaron aleatoriamente
a los 5 tratamientos que correspondían a 0, 25, 75, 200 o 500 paseos a lo largo del trayecto de una persona de
70Kg. Se trata de un diseño balanceado con 5 réplicas por cada nivel.
Los datos de este problema 1 son:
y<- c(20.7, 15.9, 17.8, 17.6, 12.9, 13.4, 12.7, 9.0,11.8, 12.6, 11.4, 12.1,7.6, 9.5,
9.9, 9.0,7.8, 9.0, 8.5, 6.7)
fact<-as.factor(c(rep(0,4),rep(25,4),rep(75,4),rep(200,4),rep(500,4)))
df<-data.frame(y,fact)

En las unidades anteriores del curso nos preguntábamos si el factor Number of passes era significativo, en el
sentido que la altura de la vegetación se vea afectada por dicho factor. Planteábamos el contraste de hipótesis

H0 : α1 = α2 = · · · = α5 = 0
H1 : αi ̸= αj para algún i ̸= j

Para resolverlo calculábamos la tabla ANOVA

results<-lm(y~fact,data=df)

Extrajimos la tabla ANOVA mediante

taov<-anova(results)
taov

## Analysis of Variance Table

##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## fact 4 243.162 60.790 29.484 5.988e-07 ***
## Residuals 15 30.927 2.062

2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Obtuvimos un pvalor 5.9877787 × 10−7 , claramente inferior a un nivel de significación del 5%. Por tanto,
aceptamos la hipótesis alternativa, y concluimos que el factor es significativo.
Recordamos que las medias de cada grupo yi. son
tapply(y,fact,mean,data=df)

## 0 25 75 200 500
## 18.000 12.000 11.975 9.000 8.000
Para el cálculo de los residuos, recordamos que cada residuo se calcula mediante rij = yij − ȳi. . Por ejemplo,
el residuo de la observación correspondiente a la primera réplica del primer nivel, es decir, i = 1 y j = 1 será

y11 − ȳ1. = 20.7 − 18 = 2.7

Podemos extraer con la función residuals todos los residuos del objeto results, de tal manera que para
cada observación obtenemos su residuo
residuals(results)

## 1 2 3 4 5
## 2.700000e+00 -2.100000e+00 -2.000000e-01 -4.000000e-01 9.000000e-01
## 6 7 8 9 10
## 1.400000e+00 7.000000e-01 -3.000000e+00 -1.750000e-01 6.250000e-01
## 11 12 13 14 15
## -5.750000e-01 1.250000e-01 -1.400000e+00 5.000000e-01 9.000000e-01
## 16 17 18 19 20
## -1.110223e-16 -2.000000e-01 1.000000e+00 5.000000e-01 -1.300000e+00

Evaluando la normalidad
Un qqplot es un procedimiento gráfico para la evaluación de la normalidad. Los qqplot son útiles en el
diagnóstico de desviaciones de la normalidad porque permiten la comparación de la cuantiles de la distribución
de los residuos de los datos con los cuantiles de una distribución normal estándar. Si se aprecian desviaciones
de la linealidad ello será indicio que la normalidad no se cumple.
Podemos obtener un qqplot de los residuos mediante la función qqnorm
qqnorm(residuals(results))
qqline(residuals(results))

3
 Normal Q−Q Plot

2
Sample Quantiles

1
0
−1
−2
−3

−2 −1 0 1 2

Theoretical Quantiles

Observamos que la mayoría de residuos se ajustan a la recta, por tanto no hay evidencia en contra del supuesto
de normalidad. Alternativamente, podemos obtener el qqplot utilizando la función plot.lm, llamando al slot
2, es decir, plot(results,which=2).
Adicionalmente podemos contrastar de normalidad mediante el test de Shapiro-Wilk. En este contraste la
hipótesis nula afirma que la distribución es normal, la hipótesis alternativa es su negación. La sentencia para
calcular este test con R es:
shapiro.test(residuals(results))

##
## Shapiro-Wilk normality test
##
## data: residuals(results)
## W = 0.96243, p-value = 0.5934
A la vista del pvalor obtenido, mantenemos la hipótesis nula y podemos aceptar que la variable aleatoria eij
sigue una distribución normal. Recordad que debemos aplicar la función shapiro.test a los residuos y es
un error (desgraciadamente muy habitual) aplicarlo sobre las observaciones yij .

Evaluando el cumplimiento de la constancia de la varianza (homocedasticidad)

Evaluamos la constancia de la varianza mediante un plot de los residuos rij en el eje vertical sobre los valores
ajustados (fitted values) en el eje horizontal. El valor ajustado de una observación es, en el anova de un
factor, la media del grupo al que corresponde la observación. Por tanto, para todos las observaciones que
corresponden nivel i, el valor ajustado es ȳi. . Por ejemplo, con los datos del ejemplo, todas las observaciones
que corresponden al nivel 2, el valor ajustado es la media del nivel 2, ȳ2. = 12.
En este tipo de gráfico de residuos aparecerán alineaciones verticales de puntos, concretamente, una tira por
cada nivel. Si la varianza es constante la extensión vertical de las tiras será aproximadamente la misma. Si
la varianza no es constante, en algunos casos aparecerá un patrón en la dispersión vertical de los residuos.
En general, en presencia de heterocedasticidad los niveles de tratamiento mostrarán una disposición de los
residuos con distinta dispersión vertical.

4
Un patrón habitual cuando no se cumple la suposición de varianza constante es aquel que aparece en caso
que la varianza depende de la media del grupo. Por lo general, en estos casos hay un aumento de la varianza
a medida que la media del grupo aumenta. En consecuencia, el gráfico de residuos muestra una apertura a la
derecha en forma de embudo (ver figura 6.3 libro).
Sigamos analizando los datos del ejemplo. Podemos obtener un gráfico de residuos utilizando la función
plot.lm, llamando al slot 1.
plot(results,which=1)

Residuals vs Fitted
3

1
2
1
Residuals

0
−1
−2

2
−3

8 10 12 14 16 18

Fitted values
lm(y ~ fact)

Observamos 4 tiras verticales de puntos, que están situadas en las medias de cada grupo, como hemos dicho
éstas corresponden a los valores ajustados de las observaciones. La disposición de los residuos muestra una
dispersión parecida en cada tira, quizá las tiras correspondientes a los niveles 500 y 200 con medias 8 y 9
respectivamente, tienen, en comparación, menor dispersión que en los otros niveles. Se aprecia por tanto un
leve efecto de embudo, pero muy débil.
Adicionalmente, podemos testar la homogeneidad de varianzas mediante el test de Levene. La hipótesis nula
es la de homogeneidad de varianzas a través de los niveles. Esto es
H0 : σ12 = σ22 = · · · = σa2
H1 : alguna diferencia

Sigamos analizando los datos del ejemplo. Evaluamos la homogeneidad de las varianzas.
library(car)

## Loading required package: carData

leveneTest(y~fact,data=df)

## Levene's Test for Homogeneity of Variance (center = median)

## Df F value Pr(>F)
## group 4 0.4262 0.7874
## 15
A la vista del pvalor, mantenemos la hipótesis nula y aceptamos que las varianzas son iguales.

5
Corrigiendo el incumplimiento de las premisas
Corrigiendo la ausencia de normalidad
La ausencia de normalidad, cuando es debida a la asimetría, puede a veces ser disminuida mediante la
transformación de la variable respuesta a una escala diferente. La asimetría hacia la derecha se ve reducida
por una raíz cuadrada, un logaritmo, u otra transformación a una potencia menor que uno, mientras que la
asimetría a la izquierda se ve reducida por un cuadrado, cubo, u otra transformación a una potencia mayor
que uno.

Corrigiendo la no constancia de la varianza

La forma habitual de solucionar la no constancia de la varianza es mediante la transformación de la variable
respuesta. Para algunas distribuciones hay transformaciones estándar que igualan o estabilizan la varianza.
Hay una teoría general de las transformaciones estabilizadoras de la varianza que se aplica a las distribuciones
donde la varianza depende de la media (Tabla 6.3, libro de Oehlert).
Cuando no disponemos una distribución con una transformación conocida para la estabilización de la varianza
(y por lo general no la tendremos) la forma tradicional de proceder ensaya una familia de transformaciones
de potencia. El método de Box-Cox establece un procedimiento para determinar a partir de los datos la
transformación de potencia. La familia de transformaciones de Box-Cox vienen dadas por
( λ
y −1
(λ) λẏ λ−1
λ ̸= 0
y =
ẏ log(y) λ = 0

donde ẏ es la media geométrica de la variable respuesta.

La técnica de Box-Cox transforma los datos (la variable respuesta) explorando un rango de valores de λ y
realizando un ANOVA para cada transformación (es decir, para cada valor de λ). En cada ANOVA obtenemos
una suma de cuadrados del error (SSE) que dependerá de la λ usada, SSE(λ). La mejor transformación λ∗
será aquella para la cual el SSE(λ) alcanza el mínimo.
Como hemos concluido anteriormente para los datos del ejemplo hemos aceptado la constancia de la varianza.
No obstante, solamente para ilustrar la aplicación de la transformación de Box-Cox, vamos a utilizarla en
estos datos. Para implementar la técnica de Box-Cox utilizaremos la función boxcox del paquete MASS.
library(MASS)
bct <- boxcox(y~fact,lambda = seq(-1, 1, length = 10),data=df,plotit=T)

6
 14.0
log−Likelihood

13.0

95%
12.0

−1.0 −0.5 0.0 0.5 1.0

lambda <- bct$x[which.max(bct$y)]

lambda

## [1] 0.3939394
Observamos la gráfica que proporciona la función boxcox donde se representa el log de la verosimilitud (en
lugar de SSE(λ)). En esta representación deberemos elegir la λ que corresponda al máximo del log de la
verosimilitud. El valor es λ∗ = 0.3939394. A nivel práctico podemos establecer para la transformación la
potencia λ = 0.4. Que no es signifcativo visto el intervalo de confianza que incluye el cero. Por otro lado, tal
como se esperaba atendiendo a los resultado de los test de normalidad y homocedasticidad. No obstante,
veamos de qué manera procederíamos en caso de ser significativo.
Transformaremos la variable respuesta de acuerdo con una λ de 0.4. Por otro lado, utilizaremos la función
geometric.mean del paquete psych para el cálculo de la media geométrica de la respuesta.
library(psych)

##
## Attaching package: 'psych'
## The following object is masked from 'package:car':
##
## logit
gm<-geometric.mean(y)
lambda<-0.4
df$yt<-(yˆlambda-1)/(lambda*gmˆ(lambda-1))

Puede ser interesante visualizar los efectos de la transformación de manera explícita comparando los valores
antes y después de la transformación.
df[,c("y","yt")]

## y yt
## 1 20.7 25.22658
## 2 15.9 21.62983

7
## 3 17.8 23.12246
## 4 17.6 22.96999
## 5 12.9 19.03686
## 6 13.4 19.49245
## 7 12.7 18.85166
## 8 9.0 15.05043
## 9 11.8 17.99582
## 10 12.6 18.75840
## 11 11.4 17.60286
## 12 12.1 18.28532
## 13 7.6 13.36731
## 14 9.5 15.61312
## 15 9.9 16.05061
## 16 9.0 15.05043
## 17 7.8 13.61855
## 18 9.0 15.05043
## 19 8.5 14.46865
## 20 6.7 12.18462
plot(df$y,df$yt,xlab="y",ylab="y transformada")
24
22
y transformada

20
18
16
14
12

8 10 12 14 16 18 20

A continuación podemos analizar de nuevo los datos pero con la respuesta transformada y hacer un gráfico
con los nuevos residuos
results2<-lm(yt~fact,data=df)
anova(results2)

## Analysis of Variance Table

##
## Response: yt
## Df Sum Sq Mean Sq F value Pr(>F)
## fact 4 212.752 53.188 27.689 9.023e-07 ***
## Residuals 15 28.814 1.921

8
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(results2,which=1)

Residuals vs Fitted

1
2
1
Residuals

0
−1

13
−2
−3

14 16 18 20 22

Fitted values
lm(yt ~ fact)

Parece que ahora la dispersión de los residuos es más similar para todos los valores ajustados en comparación
con el gráfico de residuos que obtuvimos sin transformar la respuesta. Ahora, los residuos de los niveles 500 y
200 tienen una dispersión más parecida a la de los otros niveles, se atenúa el efecto de embudo.

9
Unidad 5. Metodologías para la selección del tamaño muestral
Un aspecto importante del diseño de un experimento es determinar cuántas observaciones se necesitan para
que las conclusiones tengan una cierta precisión y una confianza suficiente. Como es lógico esperar, el tamaño
de muestra necesario queda supeditado al tipo de experimento previsto, a la forma en que se llevará a cabo y
de los recursos disponibles, de la sensibilidad deseada y la confianza exigida.
Por sensibilidad se entiende la diferencia en las medias que el experimentador desea detectar, es decir,
el experimento debe ser suficientemente sensible para detectar diferencias importantes entre los distintos
tratamientos.
En general, el aumento del número de réplicas aumenta la sensibilidad y hace que sea más fácil de detectar
pequeñas diferencias entre las medias. Tanto la potencia como el error de estimación son función del tamaño
de la muestra y función de la varianza del error.

Potencia de un contraste
El objetivo es contrastar la hipótesis de que los efectos de los distintos tratamientos son iguales, frente a la
alternativa de que algunos no lo son. La hipótesis nula de que los efectos son todos iguales se expresa como

H0 : α1 = · · · = αa = 0
Pa
ya que imponemos la restricción i=1 αi = 0 e implica que los αi son todos iguales a 0.
Como sabemos, al realizar un contraste podemos cometer un error de tipo I con probabilidad

α = P [Rechazar H0 |H0 es cierta]

o un error de tipo II si
β = P [Aceptar H0 |H1 es cierta]
La potencia de un contraste se define como

1 − β = P [Rechazar H0 |H1 es cierta]

es, por tanto, la probabilidad complementaria al error del tipo II. Se trata de construir contrastes que tengan
un nivel de significación fijo que controle α y una potencia máxima (es decir un error tipo II pequeño).
En este marco queremos calcular la potencia del test F en el caso de efectos fijos, que son estudiados hasta el
momento:

1 − β = P [F > Fα (a − 1, N − a)|H1 es cierta]

M SA
Para calcular esta probabilidad se necesita conocer la distribución de F = M SE cuando la hipótesis nula es
falsa. Se puede demostrar que en ese caso, el estadístico F se distribuye con una distribución F no centrada
con a − 1 y N − a grados de libertad y un cierto parámetro de centralidad.
Considérese un conjunto de valores de los parámetros α1 , ..., αa con los cuales H0 no es verdadera. La
probabilidad de un error de tipo II, β, es la probabilidad de que H0 no sea rechazada cuando ese conjunto es
el conjunto de valores verdaderos. Se podría pensar que β tendría que ser determinada por separado para
2
Pa de2 αi2. Afortunadamente, como β para la prueba F depende de las αi y σ sólo
cada configuración diferente
mediante la expresión i=1 αi /σ se puede evaluar al mismo tiempo para muchas alternativas diferentes.
Pa
Por ejemplo, i=1 αi2 = 4 para cada uno de los siguientes conjuntos de αi con los cuales H0 es falsa, así que
β es idéntica para las tres alternativas:
1. α1 = −1,
√ α2 = −1,
√ α3 = 1, α4 = 1
2. α1 = −√2, α2 = p2, α3 = 0, α
p4 = 0 p
3. α1 = − 3, α2 = 1/3, α3 = 1/3, α4 = 1/3

10
La cantidad
a
X
n αi2 /σ 2
i=1

se llama parámetro de no centralidad para ANOVA unidireccional. Esto es debido a que si H0 es falsa el
estadístico de prueba tiene una distribución F no centralizada con éste como uno de sus parámetros. La
probabilidad del error tipo II, β, es una función decreciente del valor de este parámetro. Por lo tanto, con
2
Pade σ2 y n, es más probable que la hipótesis nula sea rechazada para alternativas alejadas de H0
valores fijos
(es decir i=1 αi grande) que para alternativas próximas a H0 .
Pa
Con un valor fijo de i=1 αi2 , β se reduce a medida que el tamaño de muestra n en cada tratamiento se
incrementa y aumenta a medida que la varianza σ 2 se incrementa (puesto que una variabilidad subyacente
más grande dificulta detectar cualquier alejamiento dado con respecto a H0 ).

Cálculo de la potencia
Una vez que se especifican mediante sus estimaciones los valores de σ 2 y las αi para los cuales se desea β,
éstos se utilizan para calcular el valor de ϕ2
Pa 2
2 i=1 αi
ϕ =
aσ 2
ϕ recibe el nombre de medida del efecto.
Tradicionalmente se disponía de un conjunto apropiado de curvas de potencia de referencia y se escogía en
estos gráficos el tamaño de muestra. En la actualidad podemos usar por ejemplo el paquete pwr de R para
hallar tanto la potencia del test anova como el tamaño de muestra.
Usualmente se escoge una configuración de ϕ2 que recoge la mínima diferencia entre dos medias poblacionales
que el investigador considera importante. Sin perder generalidad, supongamos que son µ1 y µ2 . El resto
de a − 2 grupos en esta configuración tienen el mismo promedio µi que la media general, es decir µi = µ
para i > 2. Si se designa a esta mínima diferencia significativa entre los dos primeros grupos como 2 × δ,
tendremos que µ1 = µ + δ y µ2 = µ − δ y la configuración de las αi será:

α1 = +δ, α2 = −δ, α3 = · · · = αa = 0

de forma que la expresión del parámetro ϕ se simplifica a:

2 δ2
ϕ2 =
a σ2
En resumen, hay cuatro parámetros a escoger para calcular el tamaño de muestra del experimento: el número
de grupos a, el nivel de significación α, la potencia del contraste 1 − β y la mínima diferencia significativa 2 δ.
Hay un cierto consenso general sobre como elegir α y 1 − β (¡qué decir del tradicional nivel de significación
del 0.05!), en cambio, en la elección de δ el investigador deberá guiarse por su conocimiento experto del
campo concreto de aplicación de su estudio.
A continuación, veamos un ejemplo en el que aplicaremos los conceptos anteriores.

Ejemplo: conchas de mejillón

Los datos hacen referencia a una medición de la concha en el mejillón Mytilus trossulus, concretamente de
la longitud del músculo aductor anterior, estandarizada dividiendo por la longitud total. Se realizan en cinco
localizaciones geográficas: Tillamook, Oregon; Newport, Oregon; Petersburgo, Alaska; Magadan, Rusia; y
Tvarminne, Finlandia Los datos obtenidos aparecen en la tabla siguiente.

11
 Localización
Tillamook 0.0571 0.0813 0.0831 0.0976 0.0817 0.0859 0.0735 0.0659
Newport 0.0873 0.0662 0.0672 0.0819 0.0749 0.0649 0.0835 0.0725
Petersburgo 0.0974 0.1352 0.0817 0.1016 0.0968 0.1064 0.105 0.0829
Magadan 0.1033 0.0915 0.0781 0.0685 0.0677 0.0697 0.0764 0.0689
Tvarminne 0.0703 0.1026 0.0956 0.0973 0.1039 0.1045 0.0956 0.0818
¿Afecta de manera significativa la localización en la longitud del musculo aductor? Tomar α = 0.05.
df<-data.frame(y=c(0.0571,0.0813,0.0831,0.0976,0.0817,0.0859,0.0735,0.0659,
0.0873, 0.0662 ,0.0672,0.0819,0.0749,0.0649, 0.0835,0.0725,
0.0974 ,0.1352, 0.0817,0.1016 ,0.0968 , 0.1064,0.105 , 0.0829,
0.1033 ,0.0915 ,0.0781,0.0685 ,0.0677,0.0697 ,0.0764 , 0.0689,
0.0703,0.1026 , 0.0956, 0.0973,0.1039 ,0.1045 ,0.0956, 0.0818),
b=as.factor(c(rep(1,8), rep(2,8), rep(3,8),rep(4,8),rep(5,8) ) ))
str(df)

## 'data.frame': 40 obs. of 2 variables:

## $ y: num 0.0571 0.0813 0.0831 0.0976 0.0817 0.0859 0.0735 0.0659 0.0873 0.0662 ...
## $ b: Factor w/ 5 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 2 2 ...
boxplot(df$y~df$b)
0.12
0.10
df$y

0.08
0.06

1 2 3 4 5

df$b

results<-aov(y~b,data=df)
plot(results, which=c(1,2))

12
 Residuals vs Fitted

18

0.03 25
Residuals

0.01
−0.01

33
−0.03

0.075 0.080 0.085 0.090 0.095 0.100

Fitted values
aov(y ~ b)

Normal Q−Q
3

18
Standardized residuals

25
2
1
0
−1
−2

33

−2 −1 0 1 2

Theoretical Quantiles
aov(y ~ b)

Podemos aceptar que se cumplen las suposiciones de normalidad y homocedasticidad. Para contrastar si el
tratamiento es significativo, resolveremos la tabla anova.
anova(results)

## Analysis of Variance Table

##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## b 4 0.0042417 0.00106043 6.4605 0.0005292 ***

13
## Residuals 35 0.0057449 0.00016414
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
A la vista del pvalor aceptamos que el tratamiento es significativo. Naturalmente, que podríamos plantear
la comparación por parejas para determinar las diferencias entre tratamientos, pero a continuación nos
centraremos en el estudio de la potencia.

14
Potencia y número de réplicas
Podemos utilizar la función pwr.anova.test() del paquete pwr.
library(pwr)

En esta función intervienen los siguientes argumentos según el help de la función:

1. k, number of groups (en nuestra notación a ≡ k)
2. n, number of observations (per group)
3. f, effect size (en nuestra notación ϕ ≡ f )
4. sig.level, significance level (Type I error probability)
5. power, power of test (1 minus Type II error probability)
Exactamente uno de estos parámetros se debe pasar como NULL y entonces ese parámetro no especificacdo se
calculará en función de los otros.
Supongamos que queremos calcular la potencia de una test con la hipótesis alternativa que dos efectos difieran
de la media general en 0.01 unidades.
El tamaño del efecto ϕ serà

0.012 + (−0.01)2 + 0 + 0 + 0
ϕ2 = = 0.244
5 · 0.00016414
Y por tanto, ϕ = 0.494.
pwr.anova.test(f=0.494,k=5,n=8,sig.level=0.05)

##
## Balanced one-way analysis of variance power calculation
##
## k = 5
## n = 8
## f = 0.494
## sig.level = 0.05
## power = 0.6371954
##
## NOTE: n is number in each group
Con lo que resulta una potencia del 64% aproximadamente.
Para conocer que tamaño muestral nos proporcionaría una potencia del 80%, haremos:
pwr.anova.test(f=0.494,k=5,power = 0.8,sig.level=0.05)

##
## Balanced one-way analysis of variance power calculation
##
## k = 5
## n = 10.77398
## f = 0.494
## sig.level = 0.05
## power = 0.8
##
## NOTE: n is number in each group
Entonces, para tener una potencia del 80% seria necesario 11 replicas en lugar de las 8 actuales.

15
Test de Kruskal-Wallis: homocedasticidad sin normalidad
En situaciones en que no se pueda asumir la normalidad pero sí la homocedasticidad algunos métodos no
paramétricos nos proporcionan una herramienta para contrastar la diferencia de medias. El contraste de
Kruskal-Wallis tiene las siguientes propiedades:
1. El Objetivo del test es comprobar si hay diferencias significativas entre a grupos experimentales.
2. Kruskal-Wallis es un test no paramétrico que no asume normalidad los datos.
3. La hipótesis nula postula que la variable medida en los diferentes grupos presenta la misma distribu-
ción.
4. Asume que las formas de las distribuciones son las mismas en los diferentes grupos. Por tanto es sensible
a situaciones con datos heterocedásticos.
5. Generaliza la prueba U de Mann-Whitney para 3 o más grupos.
Podemos resumir el procedimiento del test de Kruskal-Wallis de la manera siguiente.
Procedimiento
1. Ordenar los datos de toda la muestra de menor a mayor.
2. Asignar al menor el rango = 1, al segundo rango = 2, etc. El rango de la observación i,j lo designamos
con el símbolo gij
3. El estadístico es el siguiente:

a
12 X
K= ri (ḡi. − ḡ)
N (N + 1) i=1
4. K sigue una distribución Chi-cuadrado con a-1 grados de libertad, siempre que todo ri ≥ 5.
5. En caso de empates es necesario aplicar una corrección
En R se dispone de la función kruskal.test() que podemos aplicar al ejemplo anterior.
kruskal.test(y~b,data=df)

##
## Kruskal-Wallis rank sum test
##
## data: y by b
## Kruskal-Wallis chi-squared = 15.898, df = 4, p-value = 0.003159

16
Test de Welch generalizado: normalidad con heterocedasticidad
Una generalización del test de Welch para más de dos poblaciones.
library(WRS2)

t1way(y~b,data=df, tr = 0.0, alpha = 0.05, nboot = 999)

## Call:
## t1way(formula = y ~ b, data = df, tr = 0, alpha = 0.05, nboot = 999)
##
## Test statistic: F = 5.6586
## Degrees of freedom 1: 4
## Degrees of freedom 2: 17.27
## p-value: 0.00426
##
## Explanatory measure of effect size: 0.72
## Bootstrap CI: [0.57; 0.91]
# comparacions 2 a 2
lincon(y~b,data=df, tr = 0.0, alpha = 0.05)

## Call:
## lincon(formula = y ~ b, data = df, tr = 0, alpha = 0.05)
##
## psihat ci.lower ci.upper p.value
## 1 vs. 2 0.00346 -0.01443 0.02136 0.96925
## 1 vs. 3 -0.02261 -0.04696 0.00174 0.06312
## 1 vs. 4 0.00025 -0.02052 0.02102 0.96925
## 1 vs. 5 -0.01569 -0.03570 0.00432 0.11522
## 2 vs. 3 -0.02607 -0.04886 -0.00329 0.02573
## 2 vs. 4 -0.00321 -0.02157 0.01514 0.96925
## 2 vs. 5 -0.01915 -0.03646 -0.00184 0.02652
## 3 vs. 4 0.02286 -0.00171 0.04744 0.06312
## 3 vs. 5 0.00692 -0.01715 0.03100 0.96925
## 4 vs. 5 -0.01594 -0.03630 0.00443 0.11522

17
Test de permutaciones
En ocasiones, la suposición de normalidad del término aleatorio no es plausible. En este caso es imposible
aplicar la prueba F en el contraste de la significación de los efectos del tratamiento. En tales casos se pueden
aplicar técnicas basadas en el remuestreo.
El procedimiento básico es el siguiente:
1. Calcular el valor de F para los datos experimentales (llamémosle Fobs ).
2. Elegir B muestras, donde cada muestra es una permutación aleatoria de los datos originales.
(a) Para cada muestra, asignar las primeras observaciones n1 al primer tratamiento, las siguientes
observaciones n2 al segundo tratamiento, y así sucesivamente.
(b) Para cada remuestreo, calcular F a partir de los datos.
(c) Si F > Fobs incrementar un contador
3. Después de completar el remuestreo, calcular p dividiendo el resultado en el contador por B.
4. p representa la probabilidad de obtener una F tan grande como la que obtuvimos con nuestros datos
experimentales si la hipótesis nula fuera cierta.
5. Rechazar la hipótesis nula si p < α, siendo α el nivel de significación.

Efecto de la localización en la longitud del mejillón (continuación)

Resolver el problema mediante un test de remuestreo.
Para implementar un test de remuestreo en R podemos cargar el paquete de R wPerm y llamar la fun-
ción perm.oneway.anova() y mediante el argumento R indicaremos que el remuestreo se basa en 9999
permutaciones de las observaciones.
library(wPerm)
set.seed(1234)
perm.oneway.anova(df$y,df$b,trim = 0, ford = NULL, R = 9999)

Histogram of permutation F−values

Observed
2500
Frequency

1500
500
0

0 2 4 6 8

permutation F−value

18
##
##
## RESULTS OF PERMUTATION 0% TRIMMED ONE-WAY ANOVA
## BASED ON 9999 REPLICATIONS
##
## SUMMARY STATISTICS
##
## b n Mean SD
## 1 8 0.0782625 0.012536112
## 2 8 0.0748000 0.008597176
## 3 8 0.1008750 0.016672197
## 4 8 0.0780125 0.012944656
## 5 8 0.0939500 0.012004404
##
## HYPOTHESIS TEST
##
## Response Factor Trim Statistic Observed P.value
## y b 0 F.trim 6.460543 P < 0.001
Resulta un p-value inferior al nivel de significación (5%). Por tanto, rechazamos la igualdad de los efectos.
Hay diferencias significativas de la longitud del mejillon según la localización.

19
Temas
4.1 Suposiciones del modelo.
4.2 Evaluación de las suposiciones: normalidad, varianza constante, independencia.
4.3 Posibles correcciones en caso de no cumplimiento.
4.4 Efectos en caso de no cumplimiento.
5.1 Metodologías para la selección del tamaño muestral
5.2 Potencia y tamaño muestral en ANOVA
5.3 Test de Kruskal-Wallis y test ANOVA de permutaciones

Materiales de estudio
Temario del curso Referencia bibliográfica (libro Oehlert)
De 4.1 a 4.2 Páginas de 111 a 124
De 4.3 a 4.4 Páginas de 124 a 140
De 5.1 a 5.2 Páginas de 149 a 153
De 5.3 a 5.4 Páginas de 153 a 158

20