1) Factor Simple Completamente al azar Clase
Brandon S. Solano
2023-04-25
#Diseño de un Factor simple, con un arreglo completamente al azar
#Se va a realizar una comparacion entre variedades de papa
#Imaginación visual de las pacelas (donde deberian ubicarse)
xy =expand.grid(x= seq(0,5), y = seq(0,5))
plot(xy, pch=15, cex=6, asp=1)
#Un solo factor (comunmente es el genotipo)
genotipo = gl(n = 6, k = 6, length = 36,
labels = paste0 ('gen',1:6))
#se miden las papas
#variable respuesta
##Creación de datos "recolectados", 12 cada grupo, 120 Kg pesado promedio, #desviación
set.seed(123)
Ps= c(
rnorm(12, 1200, 100),
rnorm(12, 1500, 80),
rnorm(12, 1420, 90))
#aleatorizacion en la parcela
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#DISTRIBUCIÓN DE CADA UNO DE LOS GENOTIPOS
library(ggplot2)
ggplot(datos)+
aes(x,y, fill=genotipo)+
geom_tile()
#Diseño de arbol para este diseño
library(collapsibleTree)
set.seed(123)
genotipo = data.frame(
gen = gl(6, 6, 36, c( 'gen1','gen2','gen3','gen4','gen5','gen6')),
rend = rnorm(6, mean = 1, sd = 0.2)
)
genotipo## gen rend
## 1 gen1 0.8879049
## 2 gen1 0.9539645
## 3 gen1 1.3117417
## 4 gen1 1.0141017
## 5 gen1 1.0258575
## 6 gen1 1.3430130
## 7 gen2 0.8879049
## 8 gen2 0.9539645
## 9 gen2 1.3117417
## 10 gen2 1.0141017
## 11 gen2 1.0258575
## 12 gen2 1.3430130
## 13 gen3 0.8879049
## 14 gen3 0.9539645
## 15 gen3 1.3117417
## 16 gen3 1.0141017
## 17 gen3 1.0258575
## 18 gen3 1.3430130
## 19 gen4 0.8879049
## 20 gen4 0.9539645
## 21 gen4 1.3117417
## 22 gen4 1.0141017
## 23 gen4 1.0258575
## 24 gen4 1.3430130
## 25 gen5 0.8879049
## 26 gen5 0.9539645
## 27 gen5 1.3117417
## 28 gen5 1.0141017
## 29 gen5 1.0258575
## 30 gen5 1.3430130
## 31 gen6 0.8879049
## 32 gen6 0.9539645
## 33 gen6 1.3117417
## 34 gen6 1.0141017
## 35 gen6 1.0258575
## 36 gen6 1.3430130collapsibleTreeSummary(genotipo,
hierarchy = c('gen','rend'), collapsed = F)#ANALISIS DESCRIPTIVO
#Observación del comportamiento
ggplot(datos)+
aes(genotipo, Ps)+
geom_boxplot()
ggplot(datos)+
aes(genotipo,Ps)+
geom_violin(alpha = 0.5)+ geom_boxplot(width= 0.1)
#Analisis inferencial
#todo diseño tiene un MODELO 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 de mediad (efectos)
\[y_{ij}= \mu+ \tau_{i}+\epsilon_{ij}\\
i=1,2,3,4,5,6; j=1,2,3,4,5,6\\ y_{ij}= \text{Peso Seco i-esimo genotipo
y j-esima a repeticion}\\ \mu_i= \text{media de cada i-esimogentipo}\\
\tau_{i}= \text{efecto de cada genotipo} \\ \epsilon_{ij}=
\text{residual}\] #tabla de analisis de varianza (importante)
mod1= aov(Ps ~ genotipo, data = datos)
smod1= summary(mod1)
pv1= smod1 [[1]][1,5]
ifelse(pv1<0.05,'rechazo H_0','Norechazo H _0')## [1] "rechazo H_0"#F valor: alto, Bueno ya que, la variabilidad causada por genotipos es 14 veces mayor que las repeticiones o error. #P_valor (<0.5% rechazamos la H_0, >0.5% No se rechaza H_0) - Se rechaza la hipotesis nula, lo que sugiere que existen diferencias en al menos nuno de los tratamiento en cuanto al peso seco
##Los datos propircionan evidencia en contra de la hipotesis nula (a favor de la altrna),
#Estimación de los efectos # Se recomiendo sacar la media global de todos los dato
#media global
(mu= mean(datos$Ps))## [1] 1379.69#Media por genotipo
mu_1= tapply(datos$Ps, datos$genotipo, mean)
#efecto por genotipo
tau_1 = mu_1-mu
tau_1## gen1 gen2 gen3 gen4 gen5 gen6
## -134.97483 -185.56951 123.96002 82.81483 22.23648 91.53300boxplot(Ps~genotipo,datos)
points(1:6, mu_1, pch = 16, col='blue') #media de cada tratamiento
abline(h = mu, lty=2, col='red' )#media global
segments(1:6-.2, mu_1, 1:6-0.2, mu, col='green', lwd=2, lty=2)# Efectos
Var_1= tapply(datos$Ps,datos$genotipo, var)
table(Var_1)## Var_1
## 1696.08864908128 2776.50887549191 8180.61382351066 9120.29518398284
## 1 1 1 1
## 9955.62723131245 10670.0727075404
## 1 1#revision de supuesto#
\[H_0: \sigma^2_{g1},\sigma^2_{g2},\sigma^2_{g3},\sigma^2_{g4},\sigma^2_{g5},\sigma^2_{g6} \]
#residuos
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#Prueba de bartlett- Observar si las varianzas de los genotipos son iguales o no
bt=bartlett.test(mod1$residuals,datos$genotipo)
#p_valor=0.34=3.4 #No se rechaza H_0, como el p_valor es mayor al 5%, las varianzas estadisticamente son iguales\[H_0: \epsilon \sim N(0, \sigma^2_e)\]
#prueba de normalidad
shapiro.test(mod1$residuals)##
## Shapiro-Wilk normality test
##
## data: mod1$residuals
## W = 0.97311, p-value = 0.5164#p_valor= 0.5164=51.64
#Como el p_valor es mayor al 5%, no se rechaza la H_0 de normalidad, se consideran que los residuales siguen una distribución normal#ahora se necesita COMPARAR los tratamientos y observr cuales son los que estan generando la diferencia #COMPARACIÓN DE MEDIAS POSTERIOR AL ANALISIS DE VARIANZA
#prueba de maxima diferencia de Tukey
tt= TukeyHSD(mod1,'genotipo')
plot(tt, las=1)
abline(v=0, lyt=2, col='red', lwd=2)## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...): "lyt" is
## not a graphical parameter
#los "genotipos" que presentan diferencia entre ellos estan lejos del 0, si presentan relación estaran dentro del 0
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.7103650library(agricolae)
DN= 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 bplot(DN)
#CONCLUSION: PARA SELECCIONAR EL MEJOR TRATAMIENTO ESTADISTICAMENTE PUEDE SER 3,4,5,6, QUE SON LOS MEJORES, AGRONOMICAMENTE SE SELECCIONARA EL DE MAYOR INTERES#DISEÑO SIMPLE COMPLETAMENTE AL AZAR CON INCUMPLIMIENTO DE SUPUESTOS
genotipo = gl(n = 6, k = 6, length = 36,
labels = paste0 ('gen',1:6))
set.seed(123)
Ps= c(
rnorm(12, 1200, 120),
rnorm(12, 1500, 100),
rnorm(12, 1420, 250))
#aleatorizacion en la parcela
aleat =sample(36)
datos = data.frame(xy[aleat,], genotipo,Ps)
head (datos) ## x y genotipo Ps
## 15 2 2 gen1 1132.743
## 26 1 4 gen1 1172.379
## 31 0 5 gen1 1387.045
## 16 3 2 gen1 1208.461
## 20 1 3 gen1 1215.515
## 30 5 4 gen1 1405.808datos## x y genotipo Ps
## 15 2 2 gen1 1132.7429
## 26 1 4 gen1 1172.3787
## 31 0 5 gen1 1387.0450
## 16 3 2 gen1 1208.4610
## 20 1 3 gen1 1215.5145
## 30 5 4 gen1 1405.8078
## 6 5 0 gen2 1255.3099
## 11 4 1 gen2 1048.1927
## 8 1 1 gen2 1117.5777
## 22 3 3 gen2 1146.5206
## 27 2 4 gen2 1346.8898
## 7 0 1 gen2 1243.1777
## 33 2 5 gen3 1540.0771
## 17 4 2 gen3 1511.0683
## 35 4 5 gen3 1444.4159
## 18 5 2 gen3 1678.6913
## 23 4 3 gen3 1549.7850
## 2 1 0 gen3 1303.3383
## 4 3 0 gen4 1570.1356
## 13 0 2 gen4 1452.7209
## 5 4 0 gen4 1393.2176
## 34 3 5 gen4 1478.2025
## 3 2 0 gen4 1397.3996
## 9 2 1 gen4 1427.1109
## 21 2 3 gen5 1263.7402
## 36 5 5 gen5 998.3267
## 32 1 5 gen5 1629.4468
## 10 3 1 gen5 1458.3433
## 25 0 4 gen5 1135.4658
## 14 1 2 gen5 1733.4537
## 29 4 4 gen6 1526.6161
## 12 5 1 gen6 1346.2321
## 28 3 4 gen6 1643.7814
## 1 0 0 gen6 1639.5334
## 19 0 3 gen6 1625.3953
## 24 5 3 gen6 1592.1601ggplot(datos)+
aes(genotipo, Ps)+
geom_boxplot()
mod1b = aov(Ps ~ genotipo, datos)
smod1b= summary(mod1b)# TAMBIEN LO RECHAZA OJITO JAJAJJA
shapiro.test(mod1b$residuals)##
## Shapiro-Wilk normality test
##
## data: mod1b$residuals
## W = 0.98349, p-value = 0.8558bartlett.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#p_valor dentro del 5%, se rechaza la H_0, por ellos se incumple el supuesto, las varianzas no son iguales #analisis de varianza en caso de heterocedasticidad FSCA
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#NO PERMITE EXTRAER LOS RESIDUALEScuando se incumple normalidad e igualidad de varianza usar
analisis de varianza no parametrico para ESTE diseño, no sirve para otros \[H_0: R_1=R_2=R_3=R_4=R_5=R_6\] #mucho mejor usar kruskal o permultacional (comparacion de rangos promedios posterir a kruskal.vallis)
mod1d=kruskal.test(Ps,genotipo)
#esta se puede evaluar
library(PMCMR)## Warning: package 'PMCMR' was built under R version 4.2.3## PMCMR is superseded by PMCMRplus and will be no longer maintained. You may wish to install PMCMRplus instead.##
## Attaching package: 'PMCMR'## The following object is masked from 'package:agricolae':
##
## durbin.test#posthoc.kruskal.nemenyi.test(Ps, genotipo)
library(PMCMRplus)## Warning: package 'PMCMRplus' was built under R version 4.2.3kwAllPairsNemenyiTest(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#permutacional
library(RVAideMemoire)## Warning: package 'RVAideMemoire' was built under R version 4.2.3## *** Package RVAideMemoire v 0.9-81-2 ***perm.anova(Ps ~ genotipo, data = datos, nperm=400)##
|
| | 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
## 400 permutations
## Sum Sq Df Mean Sq F value Pr(>F)
## genotipo 627712 5 125542 5.3717 0.002494 **
## Residuals 701126 30 23371
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#p_valor= es menor al 5%, se rechaza la H_0 los genotipos son diferentes