Diseño Factorial simple en arreglo completamente al azar.
Caracteristicas: 1. Único factor 2. Sin razón para bloquear
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))
genotipo
[1] gen1 gen1 gen1 gen1 gen1 gen1 gen2 gen2 gen2 gen2 gen2 gen2 gen3 gen3 gen3 gen3
[17] gen3 gen3 gen4 gen4 gen4 gen4 gen4 gen4 gen5 gen5 gen5 gen5 gen5 gen5 gen6 gen6
[33] gen6 gen6 gen6 gen6
Levels: gen1 gen2 gen3 gen4 gen5 gen6
#VARIABLE RESPUESTA
set.seed(123)
PS = c(
rnorm(12, 1200, 100),
rnorm(12, 1500, 80),
rnorm(12, 1420, 90)
)
datos = data.frame(genotipo, PS)
head(datos)
aleat = sample(36)
datos= data.frame(xy[aleat,], genotipo, PS)
head(datos)
library(ggplot2)
ggplot(datos)+
aes(x,y, fill=genotipo)+
geom_tile()
ANALISIS DESCRIPTIVO (informativo y gráfico)
#CAJAS
ggplot(datos)+
aes(genotipo, PS)+
geom_boxplot()
#VIOLINES
ggplot(datos)+
aes(genotipo, PS)+
geom_violin()
ANALISIS INFERENCIAL \[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 DE EFECTOS \[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 del i-esimo genotipo y j-esima repetición}\\ \mu_i = \text{La media de cada i-esimo genotipo}\\ \epsilon_{ij} = \text{Residuales}\]
MODELO EN FORMA DE EFECTOS \[y_{ij} = \mu + \tau_i + \epsilon_{ij}\\ \mu = \text{Media global}\\ \tau_i= \text{Efecto de cada genotipo}\]
MODELO EN FORMA MATRICIAL \[ Y = X\beta + E \\ X = \text{matriz del diseño 36 filas y 7 columnas, 1 para la media y otra para cada genotipo (6)}\\ \beta = \text{vector de parametros (\mu; \tau_1; \tau_2; \tau_3; \tau_4; \tau_5; \tau_6)}\]
OTRA FORMA DE PLANTEAR LA HIPÓTESIS \[H_0 : \tau_1 = \tau_2 = \tau_3 = \tau_4 = \tau_5= \tau_6 \]
mod1= aov(PS ~ genotipo, data=datos)
summary(mod1)
Df Sum Sq Mean Sq F value Pr(>F)
genotipo 5 502508 100502 14.22 3.67e-07 ***
Residuals 30 211996 7067
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
smod1 = summary(mod1)
smod1[[1]][1,5]
[1] 3.673791e-07
pv1 = smod1[[1]][1,5]
ifelse(pv1< 0.05,'Rechazo Ho', 'NO rechazo la Ho' )
[1] "Rechazo Ho"
Como el valor del F value es de 14,22 se puede interpretar que la variabilidad causada por los genotipos es 14,22 veces más grande que la causada por el error.
Interpretando el p-value
FORMA 1: Teniendo en cuenta que el p valor dado por la tabla de anova para este experimento es menor del 5% se recha la hipótesis nula, lo que sugiere que existen diferencias en almenos uno de los tratamientos en cuanto al peso seco.
FORMA 2: Los datos proporcionan evidencia en contra de la hipótesis nula (a favor de la alterna). Lo que dice que los tratamientos no son iguales (sin embargo no permite sacar conslusiones).
ESTIMAR EEFCTOS
#media global
mu = mean(datos$PS)
#media por genotipo
mu_i = tapply(datos$PS, datos$genotipo, mean)
#efecto por 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='yellow')
segments(1:6 -.2, mu_i, 1:6- 0.2, mu , col='blue', lwd=2, lty=2)
text(1:6, rep(1700, 6), round(tau_i, 2))
text(1:6, rep(1000, 6), round(var_i, 2))
REVISIÓN DE SUPUESTOS \[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)\] La diferencia entre lo estimado y lo calculado es lo que se conoce como residual
Lo residuales deben de tener forma de distribución normal y una varianza muy parecida o incluso igual.
hist(mod1$residuals)
var_res = tapply(mod1$residuals, datos$genotipo, var)
#TEST DE BARTLETT
#PRUEBA ESTADÍSTICA QUE PERMITE VER SI LAS VARIANZAS SON IGUALES
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
Como el p valor para la prueba de igualdad de varianzas es de 34,82 que es mayor al 5, hay evidencia cientifica para determinar que las varianzas son diferentes entre sí. Como el p valor para la prueba de igualdad de varianzas es mayor al 5% estadísticamente se pueden considerar iguales.
#PRUEBA DE SHAPIRO
#NORMALIDAD DE RESIDUOS
shapiro.test(mod1$residuals)
Shapiro-Wilk normality test
data: mod1$residuals
W = 0.97311, p-value = 0.5164
como el p valor en la prueba de normalidad es de 51.64% (>5%) se considera que los residuales siguen una distribución normal
COMPARACIÓN DE MEDIAS POSTERIOR AL ANALISIS DE VARIANZA
#PRUEBA DE TUCKEY
#(máxima diferencia de tukey)
par(mar = c(4, 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 están tocados por la linea del medio (0) se puede interpretar que no hay diferencia entre sí. En el tt cuado hay valores que tienen un p valor (p adj) mayor al 5% se puede interpretar que no difieren (son iguales)
library(agricolae)
#PRUEBA DE DUNCAN
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
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.
plot(dt)
Cuando las letras son iguales se pueden considerar genotipos iguales.
DISEÑO 1 (INCUMPLIMIENTO DE SUPUESTOS)
#FACTOR
genotipo = gl(n = 6, k = 6, length = 36,
labels = paste0('gen', 1:6))
#VARIABLE RESPUESTA
set.seed(123)
PS = c(
rnorm(12, 1200, 120),
rnorm(12, 1500, 100),
rnorm(12, 1420, 250)
)
datos = data.frame(genotipo, PS)
head(datos)
NA
ggplot(datos) +
aes(genotipo, PS) +
geom_boxplot()
mod1b = aov(PS ~ genotipo, datos)
smod1b = summary(mod1b)
mod1b
Call:
aov(formula = PS ~ genotipo, data = datos)
Terms:
genotipo Residuals
Sum of Squares 627711.6 701125.9
Deg. of Freedom 5 30
Residual standard error: 152.8753
Estimated effects may be unbalanced
shapiro.test(mod1b$residuals)
Shapiro-Wilk normality test
data: mod1b$residuals
W = 0.98349, p-value = 0.8558
bartlett.test(mod1b$residuals , datos$genotipo)
Bartlett test of homogeneity of variances
data: mod1b$residuals and datos$genotipo
Bartlett's K-squared = 12.401, df = 5, p-value = 0.02969
Como se rechaza la hipótesis de igualdad de varianzas (bartlett), se incumple el supuesto, lo cual complica la interpretación.
ANALISIS DE VARIANZA PARA UN DISEÑO FACTORIAL SIMPLE EN ARREGLO COMPLETAMENTE AL AZAR, EN PRESENCIA DE HETEROCEASTICIDAD (varianzas desiguales).
Top de idealidad:
mod1c = oneway.test(PS ~ genotipo, datos)
mod1c
One-way analysis of means (not assuming equal variances)
data: PS and genotipo
F = 8.6764, num df = 5.000, denom df = 13.702, p-value = 0.0006918
(qué hacer) cuando se incumple la normalidad e igualdad de varianzas
ANÁLISIS DE VARIANZA NO PARAMÉTRICO PARA UN DISEÑO DE ARREGLO FACTORIAL SIMPLE COMPLETAMENTE AL AZAR
\[H_0: R_1 = R_2 = R_3 = R_4 = R_5 = R_6\]
#KRUSKAL TEST
#trabaja con los rangos de los valores y no con los valores reales de los datos.
mod1d= kruskal.test(PS, genotipo)
mod1d
Kruskal-Wallis rank sum test
data: PS and genotipo
Kruskal-Wallis chi-squared = 17.204, df = 5, p-value = 0.004128
Comparación de rangos promedios posterior a kruskal.wallis
library(PMCMRplus)
kwAllPairsNemenyiTest(PS, genotipo)
Pairwise comparisons using Tukey-Kramer-Nemenyi all-pairs test with Tukey-Dist approximation
data: PS and genotipo
gen1 gen2 gen3 gen4 gen5
gen2 0.998 - - - -
gen3 0.172 0.058 - - -
gen4 0.443 0.205 0.995 - -
gen5 0.807 0.533 0.883 0.993 -
gen6 0.046 0.012 0.995 0.904 0.587
P value adjustment method: single-step
alternative hypothesis: two.sided
library(FSA)
dunnTest(PS, genotipo)
Dunn (1964) Kruskal-Wallis multiple comparison
p-values adjusted with the Holm method.
rangos = rank(PS, ties.method = 'average')
rangos
[1] 4 7 16 8 9 19 11 2 3 6 15 10 27 25 21 35 28 13 29 22 17 24 18 20 12 1
[27] 32 23 5 36 26 14 34 33 31 30
boxplot(rangos ~ genotipo)
ANALISIS DE VARIANZA PERMUTACIONAL
library(RVAideMemoire)
#Simula repeticiones
perm.anova(PS ~ genotipo, data = datos, nperm = 999 )
|
| | 0%
|
|= | 1%
|
|= | 2%
|
|== | 3%
|
|=== | 4%
|
|==== | 5%
|
|==== | 6%
|
|===== | 7%
|
|====== | 8%
|
|======= | 9%
|
|======= | 10%
|
|======== | 11%
|
|========= | 12%
|
|========== | 13%
|
|========== | 14%
|
|=========== | 15%
|
|============ | 16%
|
|============= | 17%
|
|============= | 18%
|
|============== | 19%
|
|=============== | 20%
|
|================ | 21%
|
|================ | 22%
|
|================= | 23%
|
|================== | 24%
|
|================== | 25%
|
|=================== | 26%
|
|==================== | 27%
|
|===================== | 28%
|
|===================== | 29%
|
|====================== | 30%
|
|======================= | 31%
|
|======================== | 32%
|
|======================== | 33%
|
|========================= | 34%
|
|========================== | 35%
|
|=========================== | 36%
|
|=========================== | 37%
|
|============================ | 38%
|
|============================= | 39%
|
|============================== | 40%
|
|============================== | 41%
|
|=============================== | 42%
|
|================================ | 43%
|
|================================= | 44%
|
|================================= | 45%
|
|================================== | 46%
|
|=================================== | 47%
|
|==================================== | 48%
|
|==================================== | 49%
|
|===================================== | 50%
|
|====================================== | 51%
|
|====================================== | 52%
|
|======================================= | 53%
|
|======================================== | 54%
|
|========================================= | 55%
|
|========================================= | 56%
|
|========================================== | 57%
|
|=========================================== | 58%
|
|============================================ | 59%
|
|============================================ | 60%
|
|============================================= | 61%
|
|============================================== | 62%
|
|=============================================== | 63%
|
|=============================================== | 64%
|
|================================================ | 65%
|
|================================================= | 66%
|
|================================================== | 67%
|
|================================================== | 68%
|
|=================================================== | 69%
|
|==================================================== | 70%
|
|===================================================== | 71%
|
|===================================================== | 72%
|
|====================================================== | 73%
|
|======================================================= | 74%
|
|======================================================== | 75%
|
|======================================================== | 76%
|
|========================================================= | 77%
|
|========================================================== | 78%
|
|========================================================== | 79%
|
|=========================================================== | 80%
|
|============================================================ | 81%
|
|============================================================= | 82%
|
|============================================================= | 83%
|
|============================================================== | 84%
|
|=============================================================== | 85%
|
|================================================================ | 86%
|
|================================================================ | 87%
|
|================================================================= | 88%
|
|================================================================== | 89%
|
|=================================================================== | 90%
|
|=================================================================== | 91%
|
|==================================================================== | 92%
|
|===================================================================== | 93%
|
|====================================================================== | 94%
|
|====================================================================== | 95%
|
|======================================================================= | 96%
|
|======================================================================== | 97%
|
|========================================================================= | 98%
|
|========================================================================= | 99%
|
|==========================================================================| 100%
Permutation Analysis of Variance Table
Response: PS
999 permutations
Sum Sq Df Mean Sq F value Pr(>F)
genotipo 627712 5 125542 5.3717 0.001 ***
Residuals 701126 30 23371
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Rechaza la hipotesis nula y determina que los genotipos son diferentes