Diseño Factorial al Azar
Randomized factorial design
Carlos Jiménez-Gallardo
Estadístico
MSc Infórmatica Educativa
Universidad de La Frontera
carlos.jimenez@ufrontera.cl
Data Scientist
www.innovate.cl
cjimenez@innovate.cl
Introducción
Considere la siguiente pregunta: ¿ Cúal de los nuevos fertilizantes y en que tipo de suelo se obtiene un mejor Rendimiento suelo en que sea usado?, ha diferencia de un Diseño en bloques al azar ambas variables son igualmente importantes, entendiendo que, la existencia de una combinación de estos factores pudiese provocar diferencias significativas con otras combinaciones. En la pregunta de investigación la letra “y” es clave para detectar la interacción, ya que solicita que ambos factores ocurran al mismo tiempo.
Generación del diseño
El diseño en Factorial considera que existen al menos 2 factores o variables independientes, siendo una de estas más relevante que la otra.
el modelo a considera es
\(Y_{ij}= \mu + \tau_i * \beta_j + \epsilon\)
set.seed(123)
# Diseño
fert <- c("A", "B", "C")
bloq <- paste0("B", 1:2) # solo 2 Suelos
replicas <- 1:4 # 2 réplicas por Suelo
# expandir y mostrar diseño
disenio <- expand.grid(
Fertilizante = fert,
Suelo = bloq,
Replica = replicas
)
print(disenio)## Fertilizante Suelo Replica
## 1 A B1 1
## 2 B B1 1
## 3 C B1 1
## 4 A B2 1
## 5 B B2 1
## 6 C B2 1
## 7 A B1 2
## 8 B B1 2
## 9 C B1 2
## 10 A B2 2
## 11 B B2 2
## 12 C B2 2
## 13 A B1 3
## 14 B B1 3
## 15 C B1 3
## 16 A B2 3
## 17 B B2 3
## 18 C B2 3
## 19 A B1 4
## 20 B B1 4
## 21 C B1 4
## 22 A B2 4
## 23 B B2 4
## 24 C B2 4por tanto el diseño quedo conformado como
table(disenio$Suelo, disenio$Fertilizante)##
## A B C
## B1 4 4 4
## B2 4 4 4Simulacion de datos
#SIMULACION
# Medias "verdaderas" simulacion
mu_A <- 55; mu_B <- 75; mu_C <- 58
mu_trat <- setNames(c(mu_A, mu_B, mu_C), fert)
# Efectos de Bloque (nuisance)
ef_bloq <- setNames(c( +20, +40), bloq)
#generando datos
datos <- within(disenio, {
Rendimiento <- rnorm(nrow(disenio),
mean = mu_trat[Fertilizante] * ef_bloq[Suelo],
sd = 180)
})
head(datos)## Fertilizante Suelo Replica Rendimiento
## 1 A B1 1 999.1144
## 2 B B1 1 1458.5681
## 3 C B1 1 1440.5675
## 4 A B2 1 2212.6915
## 5 B B2 1 3023.2718
## 6 C B2 1 2628.7117revisamos si se genero bien los datos
table(datos$Suelo, datos$Fertilizante)##
## A B C
## B1 4 4 4
## B2 4 4 4entonces tenemos un diseño 3x2x4 (T x S x replicas)
Cálculo sobre el Modelo
Librerias
library(tidyverse)
library(Analitica)revisar tipo de variable. ya que para el anova, las variables independientes debn ser factores.
str(datos)## 'data.frame': 24 obs. of 4 variables:
## $ Fertilizante: Factor w/ 3 levels "A","B","C": 1 2 3 1 2 3 1 2 3 1 ...
## $ Suelo : Factor w/ 2 levels "B1","B2": 1 1 1 2 2 2 1 1 1 2 ...
## $ Replica : int 1 1 1 1 1 1 2 2 2 2 ...
## $ Rendimiento : num 999 1459 1441 2213 3023 ...
## - attr(*, "out.attrs")=List of 2
## ..$ dim : Named int [1:3] 3 2 4
## .. ..- attr(*, "names")= chr [1:3] "Fertilizante" "Suelo" "Replica"
## ..$ dimnames:List of 3
## .. ..$ Fertilizante: chr [1:3] "Fertilizante=A" "Fertilizante=B" "Fertilizante=C"
## .. ..$ Suelo : chr [1:2] "Suelo=B1" "Suelo=B2"
## .. ..$ Replica : chr [1:4] "Replica=1" "Replica=2" "Replica=3" "Replica=4"mod <- aov(Rendimiento ~ Fertilizante*Suelo, data = datos)
summary(mod)## Df Sum Sq Mean Sq F value Pr(>F)
## Fertilizante 2 1456454 728227 20.99 1.97e-05 ***
## Suelo 1 10115176 10115176 291.61 1.45e-12 ***
## Fertilizante:Suelo 2 314930 157465 4.54 0.0253 *
## Residuals 18 624364 34687
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1antes de analizar corroborar supuestos
Diagnóstico Básico
# -----------------------
# Diagnósticos básicos
# -----------------------
par(mfrow = c(2,2))
plot(mod) # QQ-plot (normalidad) y Scale-Location (homocedasticidad)
par(mfrow = c(1,1))Pruebas Formales
# Pruebas opcionales
shapiro.test(mod$residuals) # normalidad de residuos##
## Shapiro-Wilk normality test
##
## data: mod$residuals
## W = 0.96974, p-value = 0.6607# install.packages("car")library(Analitica)
htrat<-Levene.Test(Rendimiento ~ Fertilizante*Suelo, data = datos,center="mean") #
summary(htrat)##
## --- Homoscedasticity Test Summary ---
##
## Method applied : Levene (mean) - global by cells
## F Statistic : 1.224693
## Degrees of freedom : 5 (between), 18 (within)
## p-value : 0.3379808 ns
## Decision (alpha = 0.05): Homoscedastic (cells)
## ----------------------------------------
##
## >> ANOVA sobre desviaciones absolutas (descomposicion factorial)
## Metodo: ANOVA on |Y - cell_center| (mean), SS Type I (sequential)
##
## Term DF F p Signif Interpretacion
## Fertilizante 2 2.1115 0.1501 ns no evidence
## Suelo 1 1.8710 0.1882 ns no evidence
## Fertilizante:Suelo 2 0.0148 0.9854 ns no evidence
## Residuals 18 NA NA no evidence
## ----------------------------------------Revisión Anova
Dado que el supuesto de normalidad es aceptado y el supuesto de homocedásticidad también, entonces se procede al análisis del ANOVA. De manera contratia, debríamos estar generando una prueba no parámetrica.
summary(mod)## Df Sum Sq Mean Sq F value Pr(>F)
## Fertilizante 2 1456454 728227 20.99 1.97e-05 ***
## Suelo 1 10115176 10115176 291.61 1.45e-12 ***
## Fertilizante:Suelo 2 314930 157465 4.54 0.0253 *
## Residuals 18 624364 34687
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1el ANOVA, establece que existe diferencias significativas por cada tratamiento, tipo de suelo y más aún habría una interaccion entre el fertilizante y el suelo, indicando la posibilidad de que al menos una de estas combinaciones responde mucho mejor o peor que otra combinación.
por ende debe analizarse el proceso Post Hoc por tratamiento, suelo e interaccion (aunque si el objetivo es la interaccion, basta con esa)
Planeando el Post Hoc
Posthoc_planner(mod,compare = c("Fertilizante"),scope=c("factor"),equal_var = TRUE)## Metrica
## 1 N levels (g)
## 2 Comparasions (m)
## 3 alpha global (FWER)
## 4 alpha_B = alpha/m (Bonf)
## 5 FWER Bonf (indep.)
## 6 alpha* Sidak
## 7 FWER Sidak (indep.)
## 8 Suggested p-value adj
## 9 Suggestion post hoc
## tabla_core
## 1 3
## 2 3
## 3 0.0500
## 4 0.0167
## 5 4.9%
## 6 0.0170
## 7 5.0%
## 8 Sidak (indep.) (alpha* = 0.016952)
## 9 Tukey HSD or Holm (LSD/SNK if you accept major Type I)Posthoc_planner(mod,compare = c("Suelo"),scope = c("factor"),equal_var = TRUE)## Metrica
## 1 N levels (g)
## 2 Comparasions (m)
## 3 alpha global (FWER)
## 4 alpha_B = alpha/m (Bonf)
## 5 FWER Bonf (indep.)
## 6 alpha* Sidak
## 7 FWER Sidak (indep.)
## 8 Suggested p-value adj
## 9 Suggestion post hoc
## tabla_core
## 1 2
## 2 1
## 3 0.0500
## 4 0.0500
## 5 5.0%
## 6 0.0500
## 7 5.0%
## 8 -
## 9 Tukey HSD or Holm (LSD/SNK if you accept major Type I)Posthoc_planner(mod,compare = c("Fertilizante","Suelo"),scope = c("cell"))## Metrica tabla_core
## 1 N levels (g) 6
## 2 Comparasions (m) 15
## 3 alpha global (FWER) 0.0500
## 4 alpha_B = alpha/m (Bonf) 0.0033
## 5 FWER Bonf (indep.) 4.9%
## 6 alpha* Sidak 0.0034
## 7 FWER Sidak (indep.) 5.0%
## 8 Suggested p-value adj Sidak (indep.) (alpha* = 0.003414)
## 9 Suggestion post hoc Holm (step-down) o Tukey HSDPostHoc
r1<-TukeyTest(mod,comparar = c("Fertilizante"))
summary(r1)## =====================================
## Multiple Comparison Method Summary
## =====================================
## Method used: Tukey
##
## >> Group means:
## A B C
## 1699.418 2226.777 1709.120
##
## >> Order of means (from highest to lowest):
## [1] "B" "C" "A"
##
## >> Pairwise comparisons:
## Comparacion Diferencia Valor_Critico p_value Significancia
## A A - B 527.3591 237.6628 0.0001 ***
## A1 A - C 9.7022 237.6628 0.9940 ns
## B B - C 517.6569 237.6628 0.0001 ***plot(r1)
r2<-TukeyTest(mod,comparar = c("Suelo"))
summary(r2)## =====================================
## Multiple Comparison Method Summary
## =====================================
## Method used: Tukey
##
## >> Group means:
## B1 B2
## 1229.235 2527.642
##
## >> Order of means (from highest to lowest):
## [1] "B2" "B1"
##
## >> Pairwise comparisons:
## Comparacion Diferencia Valor_Critico p_value Significancia
## B1 B1 - B2 1298.408 159.7413 0 ***plot(r2)
r3<-TukeyTest(mod,comparar = c("Fertilizante","Suelo"))
summary(r3)## =====================================
## Multiple Comparison Method Summary
## =====================================
## Method used: Tukey
##
## >> Group means:
## A.B1 B.B1 C.B1 A.B2 B.B2 C.B2
## 1145.116 1416.419 1126.169 2253.720 3037.135 2292.072
##
## >> Order of means (from highest to lowest):
## [1] "B.B2" "C.B2" "A.B2" "B.B1" "A.B1" "C.B1"
##
## >> Pairwise comparisons:
## Comparacion Diferencia Valor_Critico p_value Significancia
## A.B1 A.B1 - B.B1 271.3038 418.5299 0.3495 ns
## A.B11 A.B1 - C.B1 18.9470 418.5299 1.0000 ns
## A.B12 A.B1 - A.B2 1108.6048 418.5299 0.0000 ***
## A.B13 A.B1 - B.B2 1892.0191 418.5299 0.0000 ***
## A.B14 A.B1 - C.B2 1146.9561 418.5299 0.0000 ***
## B.B1 B.B1 - C.B1 290.2508 418.5299 0.2833 ns
## B.B11 A.B2 - B.B1 837.3009 418.5299 0.0001 ***
## B.B12 B.B1 - B.B2 1620.7153 418.5299 0.0000 ***
## B.B13 B.B1 - C.B2 875.6523 418.5299 0.0000 ***
## C.B1 A.B2 - C.B1 1127.5517 418.5299 0.0000 ***
## C.B11 B.B2 - C.B1 1910.9661 418.5299 0.0000 ***
## C.B12 C.B1 - C.B2 1165.9031 418.5299 0.0000 ***
## A.B2 A.B2 - B.B2 783.4144 418.5299 0.0002 ***
## A.B21 A.B2 - C.B2 38.3514 418.5299 0.9997 ns
## B.B2 B.B2 - C.B2 745.0630 418.5299 0.0003 ***plot(r3)
library(tidyverse)
resumen<-datos %>%
group_by(Fertilizante,Suelo) %>%
summarize(Prom=mean(Rendimiento))## `summarise()` has grouped output by 'Fertilizante'. You can override using the
## `.groups` argument.ggplot(resumen,aes(x=Fertilizante,y=Prom,group=Suelo,color=Suelo))+
geom_line()+
geom_point()
Conclusión
El análisis factorial muestra que tanto el tipo de fertilizante como el tipo de suelo influyen significativamente en el rendimiento del cultivo.
La interacción Fertilizante × Suelo es significativa (p < 0.05), lo que indica que el efecto de los fertilizantes depende del tipo de suelo. Específicamente, hay una mejora sustancial en el rendimiento del fertilizante B en el suelo B2, causando un efecto de interacción.