##Diseño 1 - Factorial simple en arreglo completamente al azar.
(No hay bloqueo, no hay factor que cause diferencias, no se restringe aleteorización)
#Imaginando el arreglo en el campo del experimento
xy = expand.grid(x=seq(0,5), y=seq(0,5))
plot(xy, pch=15, cex=3, asp=1)
# Factor
genotipo = gl(n = 6, k = 6, length = 36,
labels = paste0('gen', 1:6))
# Variable respuesta
set.seed(123)
PS = c(
rnorm(12, 1200, 100),
rnorm(12, 1500, 80),
rnorm(12, 1420, 90)
)
aleat = sample(36)
datos = data.frame(xy[aleat,], genotipo, PS)
head(datos)
## x y genotipo PS
## 15 2 2 gen1 1143.952
## 26 1 4 gen1 1176.982
## 31 0 5 gen1 1355.871
## 16 3 2 gen1 1207.051
## 20 1 3 gen1 1212.929
## 30 5 4 gen1 1371.506
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.2.3
ggplot(datos)+
aes(x,y, fill=genotipo)+
geom_tile()
library(ggplot2)
ggplot(datos)+
aes(genotipo, PS)+
geom_boxplot()
library(ggplot2)
ggplot(datos)+
aes(genotipo, PS)+
geom_violin()
La parte más ancha indica donde hay mas datos o donde se concentran.
*Hipotesis:
\[H_0: \mu_{g_1}=\mu_{g_2}=\mu_{g_3}=\mu_{g_4}=\mu_{g_5}=\mu_{g_6}\\ H_a: H_0\text{ es falsa}\]
*Modelo:
\[y_{ij} = \mu_i + \epsilon_{ij}\\ i=1,2,3,4,5,6-:j=1,2,3,4,5,6\]
\[y_{ij} = \text{Peso Seco i-esim genotipo y j-esima repetición}\]
\[\mu_i = \text{La media de cada i-esimo genotipo}\] \[\epsilon_{ij} = \text{residuales}\] Modelo en forma de efecto
\[y_{ij} = \mu_i = \epsilon_{ij}\] \(\mu\) media global \(\tau_i\) efecto de cada genotipo
Modelo en forma matricial
\[y = x\beta + E\]
\(X\) matriz del diseño 36 filas y 7 columnas (1 columna represengtando la media y 6 una por genotipo)
\(\beta\) vector de parametros (\(\mu: \tau_1: \tau_2: \tau_3: \tau_4: \tau_5: \tau_6\))
Otra forma de plantear la hipotesis
\[H_0: \tau_1: \tau_2: \tau_3: \tau_4: \tau_5: \tau_6\] Todos los tratamientos son iguales.(las medias son todas iguales o todas las hipotesis son nulas)
mod1 = aov(PS~genotipo, data = datos)
smod1 = summary(mod1)
pv1 = smod1[[1]][1,5]
ifelse(pv1 < 0.05, 'Rechazo Ho', 'No rechazo Ho')
## [1] "Rechazo Ho"
Interpretando el p-value: Se rechaza la hipotesis nula, lo que sugieren que existen diferencias en almenos uno de los tratamientos en cuanto al peso seco. Los datos proporcionan evidencia en contra de la hipotesis nula (a favor de la alterna)
###Estimando efectos
#Media global
(mu = mean(datos$PS))
## [1] 1379.69
#Media por genotipo
mu_i =tapply(datos$PS, datos$genotipo, mean)
#Efecto pot genotipo
tau_i = mu_i - mu
tau_i
## gen1 gen2 gen3 gen4 gen5 gen6
## -134.97483 -185.56951 123.96002 82.81483 22.23648 91.53300
var_i = tapply(datos$PS, datos$genotipo, var)
boxplot(PS ~genotipo, datos, ylim=c(1000, 1800), las=1)
points(1:6, mu_i, pch=16, col='red')
abline(h = mu, lty=2, col='red')
segments(1:6-.2, mu_i, 1:6-.2, mu, col='blue', lwdd=2, lty=2)
## Warning in segments(1:6 - 0.2, mu_i, 1:6 - 0.2, mu, col = "blue", lwdd = 2, :
## "lwdd" is not a graphical parameter
text(1:6, rep(1700,6), round(tau_i, 2))
text(1:6, rep(1000,6), round(var_i, 2))
Revisión de supuestos
varianza genotipica y varianza residual pueden ser iguales, se puede comprobar o verificar con Prueba de “bartlett” prueba estadistica
\[H_0: \sigma^2_{g1}\sigma^2_{g2}=\sigma^2_{g3}=\sigma^2_{g4}=\sigma^2_{g5}=\sigma^2_{g6}\] \[H_0: \epsilon \sim N(0, \sigma^2_e)\]
hist(mod1$residuals)
tapply(mod1$residuals, datos$genotipo, var)
## gen1 gen2 gen3 gen4 gen5 gen6
## 9120.295 8180.614 9955.627 2776.509 10670.073 1696.089
#Igualdda de varianzas
bartlett.test(mod1$residuals, datos$genotipo)
##
## Bartlett test of homogeneity of variances
##
## data: mod1$residuals and datos$genotipo
## Bartlett's K-squared = 5.5895, df = 5, p-value = 0.3482
#Normalidad de residuos
shapiro.test(mod1$residuals)
##
## Shapiro-Wilk normality test
##
## data: mod1$residuals
## W = 0.97311, p-value = 0.5164
# Prueba de máxima diferencia de Tukey
par(mar=c(5,6,3,1))
tt = TukeyHSD(mod1, 'genotipo')
plot(tt, las=1)
abline(v=0, lty=2, col='red', lwd=2)
tt
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = PS ~ genotipo, data = datos)
##
## $genotipo
## diff lwr upr p adj
## gen2-gen1 -50.594675 -198.214230 97.02488 0.8995877
## gen3-gen1 258.934855 111.315300 406.55441 0.0001225
## gen4-gen1 217.789664 70.170110 365.40922 0.0012678
## gen5-gen1 157.211312 9.591757 304.83087 0.0316170
## gen6-gen1 226.507827 78.888272 374.12738 0.0007764
## gen3-gen2 309.529530 161.909976 457.14908 0.0000068
## gen4-gen2 268.384340 120.764785 416.00389 0.0000713
## gen5-gen2 207.805987 60.186433 355.42554 0.0022109
## gen6-gen2 277.102502 129.482947 424.72206 0.0000433
## gen4-gen3 -41.145191 -188.764745 106.47436 0.9557570
## gen5-gen3 -101.723543 -249.343098 45.89601 0.3163216
## gen6-gen3 -32.427028 -180.046583 115.19253 0.9841426
## gen5-gen4 -60.578352 -208.197907 87.04120 0.8096942
## gen6-gen4 8.718162 -138.901392 156.33772 0.9999711
## gen6-gen5 69.296515 -78.323040 216.91607 0.7103650
Si un intervalo contienen el 0 significa que entre los genotipos no hay diferrencia
#Prueba de Duncan
library(agricolae)
dt= duncan.test(mod1, 'genotipo', console = T)
##
## Study: mod1 ~ "genotipo"
##
## Duncan's new multiple range test
## for PS
##
## Mean Square Error: 7066.534
##
## genotipo, means
##
## PS std r Min Max
## gen1 1244.715 95.50024 6 1143.952 1371.506
## gen2 1194.121 90.44675 6 1073.494 1322.408
## gen3 1503.650 99.77789 6 1342.671 1642.953
## gen4 1462.505 52.69259 6 1414.574 1556.108
## gen5 1401.927 103.29604 6 1268.198 1532.843
## gen6 1471.223 41.18360 6 1393.444 1500.561
##
## Alpha: 0.05 ; DF Error: 30
##
## Critical Range
## 2 3 4 5 6
## 99.11886 104.16376 107.43409 109.76839 111.53078
##
## Means with the same letter are not significantly different.
##
## PS groups
## gen3 1503.650 a
## gen6 1471.223 a
## gen4 1462.505 a
## gen5 1401.927 a
## gen1 1244.715 b
## gen2 1194.121 b
plot(dt)