Universidad Nacional de Colombia
Facultad de Ciencias Agrarias
Diseño de experimentos
Documento: 1.033.786.262 y 1.022.994.341
PARCIAL I
Nota: Para los ejercicios pares se registra como U el último número de cedula de la estudiante Yamile Mora, que corresponde al número 1. Para los ejercicios impares se registra como U y T el ultimo y tercer número de cedula de la estudiante Melisa Garzón, que corresponde al número 2 y 3.
Ejercicio Nº 1
En un estudio conducido en ambiente controlado se tuvieron 72 macetas, cada una con una planta a la que a cierta edad se le midió el contenido de clorofila (índice de clorofila) con un sensor (SPAD). El total de macetas se correspondió con 9 tratamientos asociados a estrés hídrico. Se sabe que la varianza de las 72 observaciones es 823.
Solución
- Con esta información complete la tabla del ANOVA
AOV
Datos iniciales
trt=9 #tratamientos
n=72 #total de datos analizados
rep=n/trt;rep #repeticiones## [1] 8vartot=823 #varianza total
# Suma de cuadrados de los tratamientos = 60000Resultados de tabla ANOVA
dftrt=9-1;dftrt #df tratamiento## [1] 8dftotal=72-1;dftotal #df total## [1] 71dferror=trt*(rep-1);dferror #df error## [1] 63smerror=(dftotal*vartot)-6000;smerror #Suma de cuadrados del error## [1] 52433smtotal=6000+smerror;smtotal #Suma de cuadrados del total## [1] 58433cmtrt=6000/dftrt;cmtrt #Cuadrado medio de los tratamientos## [1] 750cmerror=round(smerror/dferror);cmerror #Cuadrado medio del error## [1] 832fcal=(cmtrt/cmerror);fcal## [1] 0.9014423p.valor=pf(fcal, dftrt,dferror, lower.tail = F);p.valor## [1] 0.5209762Del F calculado podemos concluir que la variabilidad se da por las repeticiones y no por los tratamientos.
- Asumir F tabulado de 2.8 y concluir.
Hipótesis
\[H_o=\mu_{trt_1}=\mu_{trt_2}=\mu_{trt_3}=\mu_{trt_4}=\mu_{trt_5}=\mu_{trt_6}=\mu_{trt_7}=\mu_{trt_8}=\mu_{trt_9}\] Se tiene que el \(f_{cal}\) es menor que el \(f_{tab}\) (0,901<2,8), por lo tanto esto indica que se debe aceptar la hipótesis nula (Ho), es decir que las medias de los índices de clorofila entre los tratamientos no presentan diferencias entre sí, lo cual se confirma con el valor de p = 0,501, que es mayor al valor de significancia establecido del 0,05 (p=0,501>0,05).
- ¿Vale la pena comparar las medias de tratamientos a posteriori del ANOVA?
Con respecto a los resultados anteriores, no sería necesario realizar una prueba de comparación de medias ya que se determinó que no había diferencias entre ellas, por lo cual se concluye que el tratamiento de estrés hídrico no tiene influencia sobre la variable respuesta índice de clorofila.
Nota: Sin embargo, se plantea realizar una prueba de comparación de medias solamente como ejercicio práctico y académico, con el fin de analizar cómo se comportaría esta variable en un escenario diferente, donde se parte del supuesto de que se tienen medias diferentes para cada tratamiento como se pide en el punto b, lo cual es una situación más asemejada a la que podemos encontrar realmente en campo, donde efectivamente el estrés hídrico si afecta la fisiología de la planta y muy posiblemente sus valores de clorofila.
- Generar unos datos y analizar el ANOVA.
Variable respuesta: índice de clorofila.
Factor I: Tratamiento de estrés hídrico
set.seed(1) #Fijar semilla
clorofila=c(rnorm(8,402,823),rnorm(8,412,823),rnorm(8,422,823),rnorm(8,432,823),rnorm(8,442,823),rnorm(8,452,823),rnorm(8,462,823),rnorm(8,472,823),rnorm(8,482,823)) #Indice de clorofila
tto = gl(8,9,72, labels=c("T1", "T2", "T3", "T4", "T5", "T6", "T7", "T8", "T9")) #Tratamientos de estrés hídrico
datos= data.frame(clorofila,tto)I) Analisis de varianza
mod1a=aov(clorofila~tto)
summary(mod1a)## Df Sum Sq Mean Sq F value Pr(>F)
## tto 7 2374204 339172 0.565 0.781
## Residuals 64 38418126 600283De acuerdo con el análisis se concluye que hay evidencia estadística para decir que las medias son iguales entre cada tratamiento.
II) Revisión de supuestos
- Los residuos del modelo tienen distribución normal
hist(mod1a$residuals)
prueban=shapiro.test(mod1a$residuals);prueban##
## Shapiro-Wilk normality test
##
## data: mod1a$residuals
## W = 0.98862, p-value = 0.7659ifelse(prueban$p.value<0.05,'Rechazo Ho','No rechazo Ho')## [1] "No rechazo Ho"Como era de esperarse, al generar los datos con una distribución normal, efectivamente se confirma que los datos tienen este comportamiento con shapiro.test
- Las varianzas de los tratamientos son estadísticamente iguales, y esto no se duda, ya que la varianza para los datos generados es la misma.
pruebav=bartlett.test(mod1a$residuals~datos$tto);pruebav##
## Bartlett test of homogeneity of variances
##
## data: mod1a$residuals by datos$tto
## Bartlett's K-squared = 5.7088, df = 7, p-value = 0.5741ifelse(pruebav$p.value<0.05,'Rechazo Ho','No rechazo Ho')## [1] "No rechazo Ho"Se acepta Ho puesto que p=0,5741>0,05, entonces se cumple con el supuesto de igualdad de varianzas.
- Hay independencia en los residuales
plot(mod1a$residuals)
Prueba Tukey
TukeyHSD(mod1a)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = clorofila ~ tto)
##
## $tto
## diff lwr upr p adj
## T2-T1 -67.483322 -1211.8960 1076.9293 0.9999996
## T3-T1 20.570258 -1123.8424 1164.9829 1.0000000
## T4-T1 -273.981524 -1418.3942 870.4311 0.9949558
## T5-T1 37.967176 -1106.4455 1182.3798 1.0000000
## T6-T1 -75.685408 -1220.0981 1068.7273 0.9999991
## T7-T1 419.773302 -724.6394 1564.1860 0.9429071
## T8-T1 41.028943 -1103.3837 1185.4416 1.0000000
## T3-T2 88.053580 -1056.3591 1232.4662 0.9999973
## T4-T2 -206.498203 -1350.9109 937.9145 0.9991587
## T5-T2 105.450498 -1038.9622 1249.8632 0.9999908
## T6-T2 -8.202086 -1152.6147 1136.2106 1.0000000
## T7-T2 487.256624 -657.1560 1631.6693 0.8822791
## T8-T2 108.512265 -1035.9004 1252.9249 0.9999888
## T4-T3 -294.551782 -1438.9644 849.8609 0.9921687
## T5-T3 17.396918 -1127.0157 1161.8096 1.0000000
## T6-T3 -96.255666 -1240.6683 1048.1570 0.9999951
## T7-T3 399.203044 -745.2096 1543.6157 0.9560314
## T8-T3 20.458685 -1123.9540 1164.8713 1.0000000
## T5-T4 311.948700 -832.4640 1456.3614 0.9889771
## T6-T4 198.296116 -946.1165 1342.7088 0.9993540
## T7-T4 693.754827 -450.6578 1838.1675 0.5563597
## T8-T4 315.010467 -829.4022 1459.4231 0.9883243
## T6-T5 -113.652584 -1258.0652 1030.7601 0.9999846
## T7-T5 381.806127 -762.6065 1526.2188 0.9653476
## T8-T5 3.061767 -1141.3509 1147.4744 1.0000000
## T7-T6 495.458711 -648.9540 1639.8714 0.8730156
## T8-T6 116.714351 -1027.6983 1261.1270 0.9999816
## T8-T7 -378.744359 -1523.1570 765.6683 0.9668273Como conclusión, con el análisis constatamos lo que se había planteado inicialmente. Con la prueba tukey observamos que no hay diferencias entre las medias, por lo tanto, el tratamiento de estrés hídrico que se utilice no afecta el resultado final.
Ejercicio Nº 2
Antes de hilar el algodón, éste debe ser procesado para eliminar las materias extrañas y la humedad. El limpiador de pelusas más común es el limpiador de pelusas tipo sierra de batería controlada. Aunque el limpiador de pelusas de motor de sierra (M1) es uno de los más efectivos, también es uno de los limpiadores que causa más daño a las fibras de algodón. Un investigador del algodón diseñó un estudio para comparar cuatro alternativas de limpieza de las fibras de algodón: M2, M3, M4 y M5. Los métodos M2 y M3 son mecánicos, mientras que los métodos M4 y M5 son una combinación mecánica y química. El investigador quiso tener en cuenta el impacto de los diferentes cultivadores en el proceso y, por lo tanto, obtuvo fardos de algodón de seis diferentes granjas algodoneras. Las granjas fueron consideradas como bloques en el estudio. Después de una limpieza preliminar del algodón, los seis fardos fueron mezclados a fondo, y luego fue procesada una igual cantidad de algodón por cada uno de los cinco métodos de limpieza de pelusas. Las pérdidas en peso (en kg) después de la limpieza de las fibras de algodón se dan en la siguiente tabla. Durante el procesamiento de las muestras de algodón, las mediciones de la granja 1 procesada por el limpiador M1 se perdieron.
Pérdidas de algodón
Solución
Caso I. Diseño desbalanceado
Diseño: Factorial simple con bloques al azar desbalanceado.
Variable respuesta: Perdida en kg de algodón.
Factor I: Método utilizado para limpiar pelusas
Bloque: Granjero
perdida=c(NA,5.54,7.67,7.89,9.27,6.75,3.53,4.15,1.97,4.39,13.05,11.20,9.79,8.97,13.44,10.26,7.21,8.27,6.12,9.13,8.01,3.24,6.75,4.22,9.20,8.42,6.45,5.50,7.84,7.13) #Pérdida en algodón
granjero=gl(6,5,30, labels = c('G1','G2','G3','G4','G5','G6')) #Factor
metodo=gl(5,1,30, labels = c('M1','M2','M3','M4','M5')) #Bloqueo
algodon=data.frame(perdida,granjero,metodo)ggplot(data = algodon, aes(x = metodo, y = perdida, color = metodo))+geom_boxplot()+theme_bw()+ylab('Pérdida(kg)')+xlab('Método')+facet_wrap(~granjero)
En el análisis grafico encontramos que las perdidas oscilan entre 2 a 13 kg entre cada método y que los resultados varían dependiendo de la granja de donde se obtuvo el algodón. En la granja 3 se obtuvieron mayores pérdidas en los metodos 1 y 5, mientras que en la granja 2 se encuentran los mejores resultados con el método 4.
ggplot(data = algodon, aes(x = metodo, y = perdida, color = metodo))+geom_boxplot()+theme_bw()+ylab('Pérdida (kg)')+xlab('Método')
En este caso los 5 metodos parecen tener diferencias en las medias. En el método 5 se detectan algunos valores atípicos. El tamaño de las cajas es similar para todos los metodos por lo que puede que no haya indicios de falta de homocedasticidad.
I) Modelo
\[y_{ij}=\mu+\tau_{i}+\beta_{j}+\epsilon_{ij}\\i=1\dots5 \\j=1\dots6\\\sum_{i=1}^3 \tau_i=0\\ \sum_{i=1}^4 \beta_j=0\]
- \(y_{ij}\) : Variable respuesta
- \(\mu\) : Media global
- \(\tau_i\) : Efecto del método (i)
- \(\beta_j\) : Efecto de los bloques (j)
- \(\epsilon_{ij}\) : Error del método y el bloqueo
Condición lateral:
- La suma de los \(\tau_i\) es 0
- La suma de los \(\beta_j\) es 0
II) Análisis de varianza
Para el análisis de varianza se parte del modelo planteado, del cual se obtienen los siguientes resultados:
mod2a=lm(perdida~metodo+granjero)
anova(mod2a)## Analysis of Variance Table
##
## Response: perdida
## Df Sum Sq Mean Sq F value Pr(>F)
## metodo 4 47.763 11.9407 8.6494 0.0003754 ***
## granjero 5 139.661 27.9322 20.2331 5.163e-07 ***
## Residuals 19 26.230 1.3805
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Como conclusión del análisis de varianza, al observar el p-valor del método (0.0003754 < 0.05), hay evidencia estadística para decir que las medias son diferentes entre cada método utilizado.
III) Revisión de supuestos
- Los residuos del modelo tienen distribución normal
prueban=shapiro.test(mod2a$residuals);prueban##
## Shapiro-Wilk normality test
##
## data: mod2a$residuals
## W = 0.97115, p-value = 0.5913ifelse(prueban$p.value<0.05,'Rechazo Ho','No rechazo Ho')## [1] "No rechazo Ho"- Las varianzas de los tratamientos son estadísticamente iguales
En este caso no se puede aplicar la prueba Bartlett porque los datos están incompletos, entonces se utiliza como alternativa la prueba oneway.test.
pruebav2=oneway.test(perdida ~metodo,data=algodon,var.equal=FALSE);pruebav2##
## One-way analysis of means (not assuming equal variances)
##
## data: perdida and metodo
## F = 1.4544, num df = 4.000, denom df = 11.809, p-value = 0.277ifelse(pruebav2$p.value<0.05,'Rechazo Ho','No rechazo Ho')## [1] "No rechazo Ho"- Hay independencia en los residuales
plot(mod2a$residuals)
IV) Coeficiente de variación
cv=100*tapply(algodon$perdida,algodon$metodo,sd)/tapply(algodon$perdida,algodon$metodo,mean);cv## M1 M2 M3 M4 M5
## NA 47.01110 28.68060 42.91226 33.95876En el método 1 no es posible saber cuál es el coeficiente de variación debido al dato faltante. En los demás metodos el valor es superior al 20%, lo que genera desconfianza en los datos.
V. ¿Afecta el orden de colocación de los efectos del modelo dentro del software R?
mod2b=lm(perdida~granjero+metodo)
anova(mod2b)## Analysis of Variance Table
##
## Response: perdida
## Df Sum Sq Mean Sq F value Pr(>F)
## granjero 5 138.30 27.6608 20.0365 5.57e-07 ***
## metodo 4 49.12 12.2799 8.8951 0.0003186 ***
## Residuals 19 26.23 1.3805
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Los valores que se obtienen cambiando el orden de los efectos alteran levemente la suma de los cuadrados, los cuadrados medios, el F-value y el p-valor, sin embargo, esto no afecta a los residuales obtenidos en ninguno de los dos modelos, ni cambia la conclusión de las diferencias entre cada método. Para comprobar que esto no afecta se utiliza el criterio de AIC, donde podemos observar que el valor es el mismo para los dos modelos.
AIC(mod2a,mod2b)## df AIC
## mod2a 11 101.3869
## mod2b 11 101.3869Caso II. Diseño balanceado
Datos
\[Dato~faltante~(M_1.G_1)= {{6.75+13.05+10.26+8.01+8.42}\over{5}}=9.298\]
perdida_m=c(9.298,5.54,7.67,7.89,9.27,6.75,3.53,4.15,1.97,4.39,13.05,11.20,9.79,8.97,13.44,10.26,7.21,8.27,6.12,9.13,8.01,3.24,6.75,4.22,9.20,8.42,6.45,5.50,7.84,7.13)
granjero_m=gl(6,5,30, labels = c('G1','G2','G3','G4','G5','G6'))
metodo_m=gl(5,1,30, labels = c('M1','M2','M3','M4','M5'))
algodon_m=data.frame(perdida_m,granjero_m,metodo_m)ggplot(data = algodon_m, aes(x = metodo_m, y = perdida_m, color = metodo_m)) +geom_boxplot() +theme_bw()+ylab('Pérdida (kg)')+xlab('Método')
Cuando se imputa el dato faltante con la media del Método 1 (5 granjas) aparecen datos atípicos en el Método 1.
I) Análisis de varianza
Para el análisis de varianza se parte del modelo planteado en el punto anterior, del cual se obtienen los siguientes resultados:
mod2c=lm(perdida_m~metodo_m+granjero_m)
anova(mod2c)## Analysis of Variance Table
##
## Response: perdida_m
## Df Sum Sq Mean Sq F value Pr(>F)
## metodo_m 4 51.150 12.7874 9.6411 0.0001637 ***
## granjero_m 5 139.364 27.8728 21.0147 2.45e-07 ***
## Residuals 20 26.527 1.3263
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Al igual que en el análisis de varianza desbalanceado, existe evidencia estadística para decir que los resultados difieren dependiendo del método y esto se constata con el cumplimiento de los supuestos.
- Al realizar una comparación de los dos modelos se observa que cambian los grados de libertad en los residuales pasando de 19 en el desbalanceado a 20 en el balanceado. Teniendo en cuenta que se parte del número total de datos que se tienen para el análisis de varianza, el dato faltante es el que hace la diferencia.
- El análisis balanceado se muestra más preciso que el análisis desbalanceado, ya que su p-valor es inferior (0.0001637<0.0003754), aunque hay que partir de que el dato incluido no es el real, sino un dato promedio de los resultados obtenidos con el mismo método en las otras granjas, por lo tanto, no va a afectar los resultados drásticamente.
- Los valores que se obtienen comparando los 2 modelos alteran levemente la suma de los cuadrados, los cuadrados medios, el F-value y el p-valor, tanto del método, como del granjero y los residuales, pero no lo suficiente como para alterar el resultado de la diferencia de medias entre cada método.
II) Revisión de supuestos
- Los residuos del modelo tienen distribución normal
prueban=shapiro.test(mod2c$residuals);prueban##
## Shapiro-Wilk normality test
##
## data: mod2c$residuals
## W = 0.97244, p-value = 0.6079ifelse(prueban$p.value<0.05,'Rechazo Ho','No rechazo Ho')## [1] "No rechazo Ho"- Las varianzas de los tratamientos son estadísticamente iguales
pruebav=bartlett.test(mod2c$residuals~algodon_m$metodo_m);pruebav##
## Bartlett test of homogeneity of variances
##
## data: mod2c$residuals by algodon_m$metodo_m
## Bartlett's K-squared = 4.8371, df = 4, p-value = 0.3044ifelse(pruebav$p.value<0.05,'Rechazo Ho','No rechazo Ho')## [1] "No rechazo Ho"- Hay independencia en los residuales
plot(mod2c$residuals)
II) Coeficiente de variación
cv=100*tapply(algodon_m$perdida_m,algodon_m$metodo_m,sd)/tapply(algodon_m$perdida_m,algodon_m$metodo_m,mean);cv## M1 M2 M3 M4 M5
## 23.52622 47.01110 28.68060 42.91226 33.95876En el método 1, asumiendo el valor promedio del mismo método con los demás granjeros se logra un CV del 23%, todos los metodos se encuentra por encima del 20%, por lo tanto, puede existir dispersión de los datos.
Ejercicio Nº 3
Use la función de R para generar de la distribución uniforme unos datos de carbono orgánico del suelo medida a 5 cm y 10 cm de profundidad. Suponga que la medida de la capa superior osciló entre 3.0 y 3.2+0.1 y de la capa inferior osciló entre 2 y 2.3+0.2. Use expand.grid para generar una ventana de observación de 0 a 100 m para la longitud y de 0 a 200 m para la latitud. Genere 50 datos en cada capa. Use la función sort.int de R para ordenar los datos de cada capa con la opción partial=25+2 dentro de la propia función sort.int. Una vez cree los datos realice algún diagrama de color (preferiblemente 3D) que permita visualizar las medidas de carbono en cada capa generadas por computadora. Compare si se encuentran diferencias en la media de carbono entre capas utilizando un nivel de confianza del 95%.
set.seed(911)
xy=expand.grid(longitud=seq(0,100,25),latitud=seq(0,200,length.out = 10))
co_5=sort.int(runif(50,3.0,3.3),partial = 27)
co_10=sort.int(runif(50,2,2.5),partial = 27)I) Grafico 3D
datos3=data.frame(Longitud=rep(xy$longitud,2),Latitud=rep(xy$latitud,2), Profundidad = rep(c(-5,-10), each = 50), col =c(co_5,co_10))
head(datos3)## Longitud Latitud Profundidad col
## 1 0 0.00000 -5 3.026622
## 2 25 0.00000 -5 3.049018
## 3 50 0.00000 -5 3.038104
## 4 75 0.00000 -5 3.029689
## 5 100 0.00000 -5 3.024939
## 6 0 22.22222 -5 3.030435plot_ly(x=datos3$Latitud, y=datos3$Longitud, z=datos3$Profundidad, type="scatter3d", mode="markers", color = datos3$col)Como los datos se organizaron de forma espacial y a su vez de forma organizada, encontramos que a medida que nos desplazamos en las coordenadas \(x\) y \(y\), los valores de carbono orgánico son diferentes, y que a mayor profundidad de la superficie del suelo empieza a disminuir el carbono orgánico.
II) Igualdad de medias
prueba3=t.test(co_5,co_10,alternative = 't',paired = TRUE,conf.level = 0,95);prueba3##
## Paired t-test
##
## data: co_5 and co_10
## t = -10327, df = 49, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 95
## 0 percent confidence interval:
## 0.9188533 0.9188533
## sample estimates:
## mean of the differences
## 0.9188533ifelse(prueba3$p.value<0.05,'Rechazo Ho','No rechazo Ho')## [1] "Rechazo Ho"Con esta prueba evidenciamos que se rechaza la hipótesis de que las medias sean iguales entre las dos muestras de carbono orgánico. A continuación, se realiza el cálculo del coeficiente de correlación de Pearson entre los valores de CO a las dos profundidades, obteniendo un coeficiente de 0.957659 , lo cual indica que el CO tiene un alta y positiva correlación entre los 5 y 10 cm. Esto tiene sentido dado que el CO está mayormente acumulado en parte superficial del suelo, donde se encuentra la materia orgánica descompuesta o en proceso de descomposición, luego, gracias a los procesos de movimiento y flujo de nutrientes y sólidos dentro del suelo, es apenas lógico que el CO tenga correlación positiva a profundidades cercanas a la superficie (10 cm), puesto que por procesos de lixiviación y percolación este puede moverse hacia mayores profundidades de los 5 cm (FAO, 2017).
prue6=cor.test(x=co_5,y=co_10,method = 'pearson',alternative = 't',conf.level = 0.95);prue6##
## Pearson's product-moment correlation
##
## data: co_5 and co_10
## t = 23.045, df = 48, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.9262016 0.9758752
## sample estimates:
## cor
## 0.957659Ejercicio Nº 4
Caso I. ANOVA
Diseño: Factorial completo \(3^2\) en arreglo completamente al azar.
Variable respuesta: Biomasa.
Factor I: Dosis de insecticida
Factor II: Numero de aplicaciones del insecticida
Tratamientos: 9
Repeticiones: 2 para cada tratamiento
D=expand.grid(F1=c(3.25,3.75,4.25),F2=c(4,5,6))
D=rbind(D,D)
set.seed(2020)
D=D[order(sample(1:18)),]
D$biomasa=sort.int(rnorm(18,3,0.3),partial=9)
D <- within(D, F1 <- factor(F1))
D <- within(D, F2 <- factor(F2))Solución
I) Modelo del diseño
\[y_{ijk}=\mu+\tau_i+\beta_j+(\tau\beta)_{ij}+\epsilon_{ijk}\\ i=1\dots3\\j=1\dots3\\k=1,2\\\sum_{i=1}^3 \tau_i=0\\ \sum_{i=1}^3 \beta_j=0\]
\(y_{ijk}\) : Variable respuesta (Biomasa)
\(\mu\) : Media general
\(\tau_i\) : Efecto del primer factor. Dosis de insecticida (i)
\(\beta_j\) : Efecto del segundo factor. N.º. de aplicaciones del insecticida (j).
\((\tau\beta)_{ij}\) : Interacción entre dosis (i) y N.º de aplicaciones del insecticida (j).
\(\epsilon_{ijk}\) : Error de la dosis de insecticida (i), N.º de aplicaciones del insecticida (j) y repeticiones (k).
\(k\) : Repeticiones de cada tratamiento, en este caso 2
II) Análisis de varianza
mod4a=aov(D$biomasa~D$F1+D$F2+D$F1*D$F2)
summary(mod4a)## Df Sum Sq Mean Sq F value Pr(>F)
## D$F1 2 0.4161 0.20804 1.979 0.194
## D$F2 2 0.3426 0.17128 1.630 0.249
## D$F1:D$F2 4 0.1635 0.04087 0.389 0.812
## Residuals 9 0.9459 0.10510Analizando el p-valor de la interacción entre los dos factores existe evidencia estadística para concluir que no hay interacción entre los factores de dosis y el número de aplicaciones del insecticida (0.812 >0.05). Que a su vez hay evidencia estadística para decir que las medias no difieren en cada factor.
1) Prueba de comparación de medias
TukeyHSD(mod4a)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = D$biomasa ~ D$F1 + D$F2 + D$F1 * D$F2)
##
## $`D$F1`
## diff lwr upr p adj
## 3.75-3.25 0.27263391 -0.2499630 0.7952308 0.3549307
## 4.25-3.25 0.35602957 -0.1665674 0.8786265 0.1934891
## 4.25-3.75 0.08339567 -0.4392013 0.6059926 0.8975539
##
## $`D$F2`
## diff lwr upr p adj
## 5-4 -0.24371366 -0.7663106 0.2788833 0.4289015
## 6-4 -0.32456967 -0.8471666 0.1980273 0.2454279
## 6-5 -0.08085601 -0.6034529 0.4417409 0.9033230
##
## $`D$F1:D$F2`
## diff lwr upr p adj
## 3.75:4-3.25:4 0.11489233 -1.1676566 1.3974413 0.9999813
## 4.25:4-3.25:4 0.37881648 -0.9037325 1.6613655 0.9452628
## 3.25:5-3.25:4 -0.33070147 -1.6132505 0.9518475 0.9734324
## 3.75:5-3.25:4 0.16365519 -1.1188938 1.4462042 0.9997380
## 4.25:5-3.25:4 -0.07038587 -1.3529349 1.2121631 0.9999996
## 3.25:6-3.25:4 -0.37253651 -1.6550855 0.9100125 0.9497229
## 3.75:6-3.25:4 -0.16388379 -1.4464328 1.1186652 0.9997353
## 4.25:6-3.25:4 0.05642012 -1.2261289 1.3389691 0.9999999
## 4.25:4-3.75:4 0.26392415 -1.0186248 1.5464731 0.9930776
## 3.25:5-3.75:4 -0.44559380 -1.7281428 0.8369552 0.8826430
## 3.75:5-3.75:4 0.04876285 -1.2337861 1.3313118 1.0000000
## 4.25:5-3.75:4 -0.18527820 -1.4678272 1.0972708 0.9993609
## 3.25:6-3.75:4 -0.48742885 -1.7699778 0.7951201 0.8301399
## 3.75:6-3.75:4 -0.27877612 -1.5613251 1.0037729 0.9902626
## 4.25:6-3.75:4 -0.05847221 -1.3410212 1.2240768 0.9999999
## 3.25:5-4.25:4 -0.70951795 -1.9920669 0.5730310 0.4788847
## 3.75:5-4.25:4 -0.21516130 -1.4977103 1.0673877 0.9981920
## 4.25:5-4.25:4 -0.44920235 -1.7317513 0.8333466 0.8784883
## 3.25:6-4.25:4 -0.75135300 -2.0339020 0.5311960 0.4171999
## 3.75:6-4.25:4 -0.54270027 -1.8252493 0.7398487 0.7487731
## 4.25:6-4.25:4 -0.32239636 -1.6049453 0.9601526 0.9769758
## 3.75:5-3.25:5 0.49435666 -0.7881923 1.7769056 0.8206001
## 4.25:5-3.25:5 0.26031560 -1.0222334 1.5428646 0.9936553
## 3.25:6-3.25:5 -0.04183504 -1.3243840 1.2407139 1.0000000
## 3.75:6-3.25:5 0.16681768 -1.1157313 1.4493667 0.9996989
## 4.25:6-3.25:5 0.38712159 -0.8954274 1.6696706 0.9389879
## 4.25:5-3.75:5 -0.23404106 -1.5165900 1.0485079 0.9968189
## 3.25:6-3.75:5 -0.53619170 -1.8187407 0.7463573 0.7588975
## 3.75:6-3.75:5 -0.32753897 -1.6100880 0.9550100 0.9748241
## 4.25:6-3.75:5 -0.10723507 -1.3897840 1.1753139 0.9999890
## 3.25:6-4.25:5 -0.30215064 -1.5846996 0.9803983 0.9841808
## 3.75:6-4.25:5 -0.09349792 -1.3760469 1.1890511 0.9999962
## 4.25:6-4.25:5 0.12680599 -1.1557430 1.4093550 0.9999606
## 3.75:6-3.25:6 0.20865273 -1.0738963 1.4912017 0.9985347
## 4.25:6-3.25:6 0.42895663 -0.8535923 1.7115056 0.9008106
## 4.25:6-3.75:6 0.22030391 -1.0622451 1.5028529 0.9978778Con esta prueba observamos que el p-valor en todos los casos es mayor al 0.05, como era de esperar, no se encuentra diferencia significativa entre ningún par de medias.
2) Gráfico de interacción
ggplot(data = D, aes(x = F2, y = biomasa, color = F2))+geom_boxplot()+theme_bw()+ylab('Biomasa')+xlab('Dosis y Aplicación')+facet_wrap(~F1)
D %>%
group_by(F1, F2) %>%
summarise(Biom_media = mean(biomasa)) -> tips3
tips3 %>%
ggplot() +
aes(x = F1, y = Biom_media, col = F2) +
geom_line(aes(group =F2)) +
geom_point()+ylab('Biomasa')+xlab('Dosis')
A través de los gráficos observamos una interacción entre las aplicaciones 4, 5 y 6 para las dosis de 3.75 y 4.25, sin embargo, para corroborar que puede ser causada por lo cercano de los valores de biomasa, realizamos una comparación de las medias por cada factor.
Con el análisis de los gráficos podemos concluir que a mayor número de aplicaciones menor será la cantidad de biomasa presente, y que a mayor dosis del insecticida disminuye la perdida de biomasa y que en la 5 aplicación de insecticida se presenta un patrón diferente al que se evalúa en las otras dos aplicaciones. De este grafico podemos concluir que el mejor tratamiento que se podría utilizar seria dosis de 4.25 con unas 4 aplicaciones en total del insecticida.
tabla=tapply(D$biomasa,list(D$F1,D$F2), mean)
addmargins(tabla,FUN='mean',quiet=TRUE)## 4 5 6 "mean"
## 3.25 2.965294 2.634592 2.592757 2.730881
## 3.75 3.080186 3.128949 2.801410 3.003515
## 4.25 3.344110 2.894908 3.021714 3.086911
## "mean" 3.129864 2.886150 2.805294 2.940436Del análisis de las medias se observa que si comparamos las medias por dosis el mejor resultado sería el de 4.25 con una media de 3.086911, y si lo analizamos por aplicación el mejor resultado es con 4 aplicaciones con una media de 3.129864. Por lo tanto, la selección del método y la aplicación seria la misma que la obtenida en el análisis gráfico. A continuación, analizamos los supuestos.
III) Revisión de supuestos
- Los residuos del modelo tienen distribución normal
prueban=shapiro.test(mod4a$residuals);prueban##
## Shapiro-Wilk normality test
##
## data: mod4a$residuals
## W = 0.99029, p-value = 0.999ifelse(prueban$p.value<0.05,'Rechazo Ho','No rechazo Ho')## [1] "No rechazo Ho"- Las varianzas de los tratamientos son estadísticamente iguales
D$trt=interaction(D$F1,D$F2)
pruebav=bartlett.test(mod4a$residuals,D$trt);pruebav##
## Bartlett test of homogeneity of variances
##
## data: mod4a$residuals and D$trt
## Bartlett's K-squared = 7.5942, df = 8, p-value = 0.4741ifelse(pruebav$p.value<0.05,'Rechazo Ho','No rechazo Ho')## [1] "No rechazo Ho"- Hay independencia en los residuales
plot(mod4a$residuals)
Se cumple con los 3 supuestos del análisis de varianza, por lo tanto, se confía en los resultados obtenidos del ANOVA.
Caso II. ANCOVA
set.seed(2020)
D$arcilla=sort(runif(18,0.20,0.40),decreasing = TRUE)I) Modelo
\[y_{ijk}=\mu+\tau_i+\beta_j+\alpha(x_{ijk}- \bar x)+(\tau\beta)_{ij}+\epsilon_{ijk}\\ i=1\dots3\\j=1\dots3\\k=1,2\\y~sus~condiciones~laterales\]
\(y_{ijk}\) : Variable respuesta (Biomasa)
\(\mu\) : Media general
\(\tau_i\) : Efecto del primer factor. Dosis de insecticida (i)
\(\beta_j\) : Efecto del segundo factor. N.º. de aplicaciones del insecticida (j).
\(\alpha\) : Es el coeficiente de regresión lineal que representa la cantidad en que varía la variable respuesta por cada cambio unitario de covariable
\(x_{ijk}\) : Valor de la covariable correspondiente a la observación \(y_{ijk}\).
\(\bar x\) : Es la media de la covariable
\((\tau\beta)_{ij}\) : Interacción entre dosis (i) y N.º de aplicaciones del insecticida (j).
\(\epsilon_{ijk}\) : Error de la dosis de insecticida (i), N.º de aplicaciones del insecticida (j) y repeticiones (k).
\(k\) : Repeticiones de cada tratamiento, en este caso 2
II) Revisión de supuestos
Para poder analizar un ANCOVA primero revisamos los supuestos
1) Linealidad entre la covariable y la variable de resultado en cada nivel de la variable de agrupación. Esto se puede verificar creando un diagrama de dispersión agrupado de la covariable y la variable de resultado.
2) Homogeneidad de las pendientes de regresión. Las pendientes de las líneas de regresión, formadas por la covariable y la variable de resultado, deben ser las mismas para cada grupo. Esta suposición evalúa que no hay interacción entre el resultado y la covariable. Las líneas de regresión trazadas por grupos deben ser paralelas.
3) La variable de resultado debe tener una distribución aproximadamente normal. Esto se puede verificar usando la prueba de normalidad de Shapiro-Wilk en los residuos del modelo.
4) Homocedasticidad u homogeneidad de la varianza de residuos para todos los grupos. Se supone que los residuos tienen una varianza constante (homocedasticidad). 5) Independencia
Tomado de: [Datanovia]
1) Linealidad
ggscatter(D, x = "arcilla", y = "biomasa", add = "reg.line",xlab = "Arcilla",ylab = "Biomasa",cor.coef = TRUE,color = "orange")
lin=cor.test(D$biomasa, D$arcilla);lin##
## Pearson's product-moment correlation
##
## data: D$biomasa and D$arcilla
## t = -2.9546, df = 16, p-value = 0.00932
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.8306127 -0.1761319
## sample estimates:
## cor
## -0.594145A través del grafico comparamos la variable respuesta con la covariable de arcilla, observamos que a medida que se aumenta el contenido de arcilla los valores de biomasa disminuyen. Para comprobar si existe una correlación se utiliza la prueba de correlación de Pearson, donde concluimos que hay evidencia estadística para decir que existe linealidad entre la covariable y la variable respuesta.
2) Homogeneidad pendientes de regresión
anova3 = aov(D$biomasa~D$arcilla*D$F1)
summary(anova3)## Df Sum Sq Mean Sq F value Pr(>F)
## D$arcilla 1 0.6594 0.6594 9.816 0.00864 **
## D$F1 2 0.2405 0.1203 1.790 0.20879
## D$arcilla:D$F1 2 0.1619 0.0810 1.205 0.33344
## Residuals 12 0.8062 0.0672
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova4 = aov(D$biomasa~D$arcilla*D$F2)
summary(anova4)## Df Sum Sq Mean Sq F value Pr(>F)
## D$arcilla 1 0.6594 0.6594 10.480 0.00712 **
## D$F2 2 0.4530 0.2265 3.599 0.05962 .
## D$arcilla:D$F2 2 0.0006 0.0003 0.004 0.99557
## Residuals 12 0.7551 0.0629
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Hay evidencia estadística para decir que hay homogeneidad de pendientes de regresión.
3) Distribución normal
anova3 = lm(D$biomasa~D$arcilla+D$F1*D$F2)
anova(anova3)## Analysis of Variance Table
##
## Response: D$biomasa
## Df Sum Sq Mean Sq F value Pr(>F)
## D$arcilla 1 0.65945 0.65945 11.6520 0.009176 **
## D$F1 2 0.24050 0.12025 2.1248 0.181922
## D$F2 2 0.46069 0.23034 4.0700 0.060359 .
## D$F1:D$F2 4 0.05468 0.01367 0.2415 0.906943
## Residuals 8 0.45276 0.05660
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1shapiro_test(anova3$residuals)## # A tibble: 1 x 3
## variable statistic p.value
## <chr> <dbl> <dbl>
## 1 anova3$residuals 0.965 0.694Hay evidencia estadística para decir que hay distribución normal
4) Independencia
anova1 = aov(D$arcilla~D$F1) #anova covariable~factor 1. Dosis
summary(anova1)## Df Sum Sq Mean Sq F value Pr(>F)
## D$F1 2 0.01603 0.008015 3.486 0.0571 .
## Residuals 15 0.03449 0.002299
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova2 = aov(D$arcilla~D$F2) #anova covariable~factor 2. Aplicaciones
summary(anova2)## Df Sum Sq Mean Sq F value Pr(>F)
## D$F2 2 0.00150 0.000748 0.229 0.798
## Residuals 15 0.04902 0.003268Hay evidencia estadística para decir que hay homogeneidad en las pendientes de regresión ya que no existe interacción entre la covariable de arcilla y los factores de dosis y numero de aplicaciones.
5) Homogeneidad de la varianza
leveneTest(D$biomasa, D$F1, center="median") #Prueba Homocedasticidad factor Dosis## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 2 0.3417 0.716
## 15leveneTest(D$biomasa, D$F2, center="median") #Prueba Homocedasticidad factor aplicacion## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 2 0.066 0.9364
## 15Hay evidencia estadística para decir que hay homogeneidad de varianzas.
III) Análisis de covarianza
ANCOVAt = aov(biomasa~arcilla+F1+F2+F1*F2, data=D) #efecto de los factores, covariable e interaccion de los factores
summary(ANCOVAt)## Df Sum Sq Mean Sq F value Pr(>F)
## arcilla 1 0.6594 0.6594 11.652 0.00918 **
## F1 2 0.2405 0.1203 2.125 0.18192
## F2 2 0.4607 0.2303 4.070 0.06036 .
## F1:F2 4 0.0547 0.0137 0.242 0.90694
## Residuals 8 0.4528 0.0566
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ANCOVAt = aov(biomasa~arcilla+F1+F2, data=D) #efecto de los factores y la covariable
summary(ANCOVAt)## Df Sum Sq Mean Sq F value Pr(>F)
## arcilla 1 0.6594 0.6594 15.595 0.00193 **
## F1 2 0.2405 0.1203 2.844 0.09752 .
## F2 2 0.4607 0.2303 5.447 0.02074 *
## Residuals 12 0.5074 0.0423
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Analizando el p-valor de la interacción entre los dos factores existe evidencia estadística para concluir que no hay interacción entre los factores de dosis y el número de aplicaciones del insecticida (0.90694 >0.05). Como la interacción no influye la retiramos del modelo y obtenemos un nuevo análisis de covarianza, donde observamos que existe evidencia estadística para rechazar la hipótesis nula de medias iguales en cada factor y covariable. En este caso la covariable influye en el contenido de biomasa al igual que los factores.
Para el ANOVA y el ANCOVA se cumplen los supuestos, por lo tanto, los resultados de los análisis son confiables.
Comparando los gráficos de interacción del ANOVA y el ANCOVA, nos damos cuenta de que el contenido de arcillas influye en el contenido de biomasa, lo cual influye en la determinación de que dosis utilizar, si se tratara de mejor biomasa se optaría por la dosis de 4.25 pero con suelos bajos en arcillas, pero también se observan resultados similares si realizo el tratamiento con una dosis de 3.75 en suelos con mayor proporción de arcillas.
Si recomendásemos el uso de arcillas expandibles, esta covariable muestra mejores resultados con dosis más bajas cuando los contenidos son altos, un ejemplo en la dosis de 3.75 que se asemeja a la dosis de 4.25 en cuanto a contenido de biomasa, esto ayudaría al agricultor a reducir costos en insecticidas si tiene un suelo con contenido de arcillas alto.
Es un hecho que el uso de insecticidas disminuye la cantidad de biomasa presente en el sustrato suelo, así que es recomendable hacer uso de menores frecuencias y dosis.
Ejercicio Nº 5
Existe un tipo de diseño anidado (factorial incompleto) conocido como anidado escalonado (staggered nested design) y ocurre tal como se muestra en la imagen, donde se tienen dos fincas sembradas con variedades de papa, solo que la finca A permite que se desarrollen las dos variedades mientras que la altitud de la finca B solo permite el desarrollo de una de ellas. Además, se tienen dos parcelas con la variedad 1 en la primera finca y solo una en el resto de las fincas.
Esquema de arbol
- Use la librería lme4 tal como aparece en el código abajo. La etiqueta “ue” hace referencia a la unidad experimental (parcela) utilizada, por lo que se necesita crear una columna que identifique la parcela, una que identifique la finca, otra para la variedad y otra para lo que aquí se llama test pero que hace referencia en este caso a los cuadrados de 1.5m*1.5m usados para tomar las muestras de plantas dentro de las parcelas. Estos diseños son usados para estimar la varianza atribuible a las parcelas, a las parcelas anidadas en las fincas, y a la variedad dentro de la finca. El código presentado puede ayudar a la estimación de estas varianzas.
Modelo
Solución
Datos
- Estimación de las varianzas
Datos
I) Modelo diseño anidado escalonado
mod5b= lmer(Respuesta ~ 1 + (1|UE) + (1|Finca) + (1|Variedad)+ (1|Repetición) + (1|UE:Finca)+
+ (1|UE:Finca:Variedad), data = data5) #modelo
summary(mod5b)## Linear mixed model fit by REML ['lmerMod']
## Formula: Respuesta ~ 1 + (1 | UE) + (1 | Finca) + (1 | Variedad) + (1 |
## Repetición) + (1 | UE:Finca) + +(1 | UE:Finca:Variedad)
## Data: data5
##
## REML criterion at convergence: 323.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.69395 -0.39641 -0.05901 0.46924 1.59881
##
## Random effects:
## Groups Name Variance Std.Dev.
## UE:Finca:Variedad (Intercept) 1.2075 1.0988
## UE:Finca (Intercept) 0.0000 0.0000
## UE (Intercept) 7.0542 2.6560
## Repetición (Intercept) 0.1492 0.3863
## Variedad (Intercept) 0.0000 0.0000
## Finca (Intercept) 0.0000 0.0000
## Residual 0.7892 0.8884
## Number of obs: 80, groups:
## UE:Finca:Variedad, 60; UE:Finca, 40; UE, 20; Repetición, 2; Variedad, 2; Finca, 2
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 8.0873 0.6807 11.88
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingularrep=(100*(0.1492)/(0.1492+7.0542+1.2075));rep## [1] 1.773889UE=(100*(7.0542)/(0.1492+7.0542+1.2075));UE## [1] 83.86974UE_fin_var=(100*(1.2075)/(0.1492+7.0542+1.2075));UE_fin_var## [1] 14.35637De acuerdo con los resultados mostrados en la tabla anterior, se tiene que la mayor proporción de la varianza se asocia a la unidad experimental (UE), es decir a las parcelas, con un valor del 83% y solo un 2% se asocia a las repeticiones, es decir a la variabilidad dentro de cada parcela, también se observa que la varianza por los factores finca y variedad es cero. Se puede concluir que este es un escenario ideal para la evaluación de diferentes clones de papa en un trabajo de fitomejoramiento, para el cual son muy usados este tipo de diseños, ya que permiten estimar la varianza asociada al ambiente para que no interfiera en el análisis cuando se evalúe el potencial de rendimiento o cualquier otra variable evaluada (Bandera y Pérez, 2018).
Para explicar mejor estos resultados, se plantea el contexto de un caso donde se esté evaluando el rendimiento de dos variedades de papa (variedad), cada una de ellas con 20 clones o genotipos distintos (UE), en dos ambientes diferentes (finca). Entonces podría plantearse que variabilidad se asocia a las condiciones de la unidad experimental evaluada en la mayor parte, es decir que cada clon de papa tiene una respuesta diferente en su rendimiento aun siendo de la misma variedad, además estos genotipos en cada variedad se comportan de manera similar en los diferentes ambientes (finca), lo que se puede suponer que son materiales estables y con buen potencial de respuesta para ser liberados en campo, previa selección de los mejores materiales en cuanto a la variable que se estaba evaluando, en este ejemplo, el rendimiento.
Ejercicio Nº 6
En el enlace [ASBIO] se tienen unos datos de potasio de muestras de suelos medidas en 8 diferentes laboratorios.
data("K")
head(K)## K lab
## 1 296 B
## 2 260 B
## 3 341 B
## 4 359 B
## 5 323 B
## 6 321 BSolución
- Análisis descriptivo
A continuación, se realiza un análisis descriptivo con medidas de tendencia central (Media, Mediana y Mediana truncada), medidas de variabilidad (Desviación estándar, Rango) y medidas de forma (curtosis) para todo el conjunto general de datos. El error estándar que aquí se menciona es un dato erróneo, ya que se calcula teniendo en cuenta una única media.
describe(K$K)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 72 307.79 48.4 314 311.1 43.74 187 413 226 -0.54 0.04 5.7Análisis Global
Media global
mean=mean(K$K);mean## [1] 307.7917Desviación estándar global
sd=sd(K$K);sd## [1] 48.39682Coeficiente de variación global
100*(sd)/(mean)## [1] 15.72389Si se realiza un análisis a nivel global, no se sospecharía de los resultados que se esconden entre e intra-laboratorio.
Análisis por laboratorio
Media por laboratorio
mean_lab=tapply(K$K,K$lab,mean);mean_lab## B D E F G H I J
## 326.1111 321.1111 316.5556 315.6667 304.1111 229.3333 313.2222 336.2222Por medio de la media observamos que hay diferencias en los resultados obtenidos por cada laboratorio. Siendo la misma muestra para todos los laboratorios los resultados oscilan entre 230 a 336 mg/kg, diferencias muy altas.
Desviación estándar por laboratorio
sd_lab=tapply(K$K,K$lab,sd);sd_lab## B D E F G H I J
## 44.40564 29.25510 80.35097 25.04995 15.37404 23.89037 31.39577 21.60890Para comparar los resultados intra-laboratorio se utiliza la desviación estándar de cada laboratorio, se asume que si las muestras son las mismas estos valores deberían ser muy bajos, por lo tanto, esto denota que los procesos que se llevan a cabo dentro del laboratorio no son controlados o parametrizados produciendo desviaciones mayores, en consecuencia, también se sigue observando diferencias entre laboratorios.
Coeficiente de variación por laboratorio
cv_lab=100*sd_lab/mean_lab;cv_lab## B D E F G H I J
## 13.616722 9.110586 25.382896 7.935570 5.055402 10.417315 10.023481 6.426969- Análisis gráfico
ggplot(data = K, aes(x = lab, y = K, color = lab)) +
geom_boxplot() +
theme_bw()+ylab('Potasio(mg/kg)')+xlab('Laboratorio')
ggplot(data=K,aes(x=K,fill=lab))+
geom_density(alpha=0.3)+ylab('Densidad')+xlab('Potasio(mg/kg)')
Analizando el diagrama de caja observamos que los tamaños de las cajas son diferentes entre cada laboratorio, por ende, se verá reflejado en un análisis de varianza, afectando la homocedasticidad y la normalidad de los residuales. También encontramos datos atípicos en los laboratorios B, D y H.
Por medio del grafico de densidades, observamos que la normalidad se alcanza en unos pocos laboratorios, la mayoría tiene distribuciones asimétricas.
A través del análisis matemático y con los gráficos, observamos que hay diferencias en las medias entre cada tratamiento y asumiendo que la muestra es la misma para todos los laboratorios, hay variabilidad de los resultados entre laboratorios y dentro de los mismos laboratorios.
Caso I. Análisis de varianza
mod5a=aov(K~lab,data=K)
summary(mod5a)## Df Sum Sq Mean Sq F value Pr(>F)
## lab 7 68930 9847 6.472 8.92e-06 ***
## Residuals 64 97370 1521
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Del análisis se obtiene que hay evidencia estadística para decir que las medias son diferentes entre cada laboratorio, a continuación, se revisan los supuestos.
- NO cumple con normalidad
prueban=shapiro.test(mod5a$residuals);prueban##
## Shapiro-Wilk normality test
##
## data: mod5a$residuals
## W = 0.94688, p-value = 0.004268ifelse(prueban$p.value<0.05,'Rechazo Ho','No rechazo Ho')## [1] "Rechazo Ho"- NO cumple con igualdad de varianzas
pruebav=bartlett.test(mod5a$residuals,K$lab);pruebav##
## Bartlett test of homogeneity of variances
##
## data: mod5a$residuals and K$lab
## Bartlett's K-squared = 32.201, df = 7, p-value = 3.727e-05ifelse(pruebav$p.value<0.05,'Rechazo Ho','No rechazo Ho')## [1] "Rechazo Ho"- Cumple con la independencia
plot(mod5a$residuals)
Si realizamos un análisis de varianza nos damos cuenta de que no se cumplen los supuestos de normalidad en los residuales ni igualdad de varianzas, por lo tanto, se recurre a la Prueba de Kruskall-Wallis, prueba no paramétrica para verificar si las medias de los laboratorios son diferentes.
Hipótesis nula
\[H_o:\mu_{rango_{lab_A}}=\mu_{rango_{lab_D}}=\mu_{rango_{lab_E}}=\mu_{rango_{lab_F}}=\mu_{rango_{lab_G}}=\mu_{rango_{lab_H}}=\mu_{rango_{lab_I}}=\mu_{rango_{lab_J}}\]
krusk=kruskal.test(K$K~K$lab);krusk##
## Kruskal-Wallis rank sum test
##
## data: K$K by K$lab
## Kruskal-Wallis chi-squared = 24.482, df = 7, p-value = 0.000937ifelse(krusk$p.value<0.05,'Rechazo Ho','No rechazo Ho')## [1] "Rechazo Ho"Con esta prueba rechazamos la hipótesis nula, lo que quiere decir que los rangos promedios asignados a cada laboratorio son diferentes. Por lo tanto, no se requiere del supuesto de normalidad e igualdad de varianzas.
Prueba para comparar por pares de laboratorio
- Prueba HSD.test y Tukey HSD
Test de comparaciones múltiples, que permite comparar los laboratorios y jerarquizarlos de acuerdo con sus medias.
TukeyHSD(mod5a)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = K ~ lab, data = K)
##
## $lab
## diff lwr upr p adj
## D-B -5.0000000 -62.61388 52.61388 0.9999939
## E-B -9.5555556 -67.16944 48.05833 0.9995151
## F-B -10.4444444 -68.05833 47.16944 0.9991328
## G-B -22.0000000 -79.61388 35.61388 0.9300514
## H-B -96.7777778 -154.39166 -39.16389 0.0000464
## I-B -12.8888889 -70.50277 44.72499 0.9966852
## J-B 10.1111111 -47.50277 67.72499 0.9992979
## E-D -4.5555556 -62.16944 53.05833 0.9999968
## F-D -5.4444444 -63.05833 52.16944 0.9999891
## G-D -17.0000000 -74.61388 40.61388 0.9825227
## H-D -91.7777778 -149.39166 -34.16389 0.0001272
## I-D -7.8888889 -65.50277 49.72499 0.9998643
## J-D 15.1111111 -42.50277 72.72499 0.9912314
## F-E -0.8888889 -58.50277 56.72499 1.0000000
## G-E -12.4444444 -70.05833 45.16944 0.9973396
## H-E -87.2222222 -144.83611 -29.60834 0.0003118
## I-E -3.3333333 -60.94722 54.28055 0.9999996
## J-E 19.6666667 -37.94722 77.28055 0.9608089
## G-F -11.5555556 -69.16944 46.05833 0.9983369
## H-F -86.3333333 -143.94722 -28.71945 0.0003704
## I-F -2.4444444 -60.05833 55.16944 1.0000000
## J-F 20.5555556 -37.05833 78.16944 0.9505152
## H-G -74.7777778 -132.39166 -17.16389 0.0031670
## I-G 9.1111111 -48.50277 66.72499 0.9996458
## J-G 32.1111111 -25.50277 89.72499 0.6574404
## I-H 83.8888889 26.27501 141.50277 0.0005920
## J-H 106.8888889 49.27501 164.50277 0.0000057
## J-I 23.0000000 -34.61388 80.61388 0.9130199HSD=HSD.test(mod5a,'lab');HSD## $statistics
## MSerror Df Mean CV MSD
## 1521.406 64 307.7917 12.6726 57.61388
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey lab 8 4.431245 0.05
##
## $means
## K std r Min Max Q25 Q50 Q75
## B 326.1111 44.40564 9 260 413 296 323 341
## D 321.1111 29.25510 9 266 354 315 326 343
## E 316.5556 80.35097 9 187 400 283 351 376
## F 315.6667 25.04995 9 274 354 305 308 327
## G 304.1111 15.37404 9 280 326 297 301 316
## H 229.3333 23.89037 9 199 280 218 226 241
## I 313.2222 31.39577 9 269 355 284 311 339
## J 336.2222 21.60890 9 299 359 318 342 353
##
## $comparison
## NULL
##
## $groups
## K groups
## J 336.2222 a
## B 326.1111 a
## D 321.1111 a
## E 316.5556 a
## F 315.6667 a
## I 313.2222 a
## G 304.1111 a
## H 229.3333 b
##
## attr(,"class")
## [1] "group"- Prueba de Pairwase t
Test de comparaciones múltiples, que permite comparar los laboratorios y hacer correcciones de ser necesario.
pairwise.t.test(K$K,K$lab)##
## Pairwise comparisons using t tests with pooled SD
##
## data: K$K and K$lab
##
## B D E F G H I
## D 1.00000 - - - - - -
## E 1.00000 1.00000 - - - - -
## F 1.00000 1.00000 1.00000 - - - -
## G 1.00000 1.00000 1.00000 1.00000 - - -
## H 4.7e-05 0.00013 0.00030 0.00035 0.00293 - -
## I 1.00000 1.00000 1.00000 1.00000 1.00000 0.00054 -
## J 1.00000 1.00000 1.00000 1.00000 1.00000 5.9e-06 1.00000
##
## P value adjustment method: holm- Prueba Nemenyi
Test de comparaciones múltiples, que permite comparar los laboratorios por rangos utilizados en la prueba no paramétrica Kruskal Wallis.
posthoc.kruskal.nemenyi.test(K~lab, data=K)## Warning in posthoc.kruskal.nemenyi.test.default(c(296, 260, 341, 359, 323, :
## Ties are present, p-values are not corrected.##
## Pairwise comparisons using Tukey and Kramer (Nemenyi) test
## with Tukey-Dist approximation for independent samples
##
## data: K by lab
##
## B D E F G H I
## D 1.0000 - - - - - -
## E 1.0000 1.0000 - - - - -
## F 0.9999 0.9999 0.9998 - - - -
## G 0.9324 0.9324 0.9222 0.9943 - - -
## H 0.0098 0.0098 0.0087 0.0397 0.2764 - -
## I 0.9993 0.9993 0.9989 1.0000 0.9984 0.0600 -
## J 0.9893 0.9893 0.9916 0.9051 0.4405 0.0003 0.8461
##
## P value adjustment method: noneDe todas las pruebas que compararon las medias y rangos se observa que el laboratorio H es el que presenta diferencias respecto a los otros laboratorios, y eso corresponde a la media más baja de la cantidad de potasio en las muestras de suelo evaluadas.
Como resumen del análisis se logra percibir que hay diferencias en los resultados de potasio determinado por cada laboratorio e inclusive el análisis dentro del laboratorio, esto ocurre usualmente cuando no se cumplen con los protocolos de calidad dentro de cada laboratorio, partiendo de que cada procedimiento tiene un protocolo para reducir al máximo errores durante el ensayo e inclusive cuando se registran los datos, se recomienda acudir a laboratorios que estén certificados en cuanto a control de calidad de sus procesos para evitar utilizar resultados inconsistentes en nuestros análisis.
Ejercicio Nº 7
Diseñe un experimento en parcelas divididas en bloques completos (diseño en franjas o strip plot design). Genere los datos usted mismo y esquematice el diseño.
- Expliqué las razones de colocar el primer factor en la parcela principal y el segundo en la subparcela.
- Genere unos datos asociados a una covariable y corra el análisis de covarianza respectivo. ¿Se justifica el uso de la covariable en el modelo? ¿Se justifica el bloque en el modelo? ¿Se tiene interacción de factores?
- De no presentarse interacción, reduzca el modelo a la presencia de solo términos cuyos p_valores sean menores al 6%. Escriba el modelo final e interprete el resultado desde un punto de vista agronómico seleccionando el “mejor tratamiento” en la mejor condición de bloqueo y con la presencia de la covariable. No olvide ordenar datos de la covariable. Revise los supuestos necesarios para el análisis estadístico que está proponiendo.
Solución
Datos
Generación de datos de rendimiento y pH: Se utiliza la función rnorm para rendimiento y runif para pH. Los parámetros de media y desviación estándar se aproximan con el criterio técnico agronómico propio según la respuesta esperada para cada tratamiento. Para el caso del pH, se definen dos intervalos para cada bloque según lo indicado anteriormente, de modo que para los datos del bloque 1 se usa un intervalo de (4.3 a 5.5) y para el bloque 2 de (3.5 a 4.3). Enseguida se crea la matriz de datos con ayuda de la herramienta excel.
I) Definición del diseño experimental
Se plantea realizar la evaluación del rendimiento de papa (Solanum tuberosum) variedad Pastusa Suprema, aplicando 8 tratamientos, a fin de determinar el efecto de cada uno de ellos y el efecto de la interacción si lo hay, dado que el cultivo de papa limita su producción cuando no se le proporciona las dosis adecuadas de agua de riego como también del fertilizante, en especial los que aportan nitrógeno. Por lo tanto, se plantean las siguientes condiciones del experimento: Se dispone de un lote, el cual se propone dividir en 2 parcelas por la disponibilidad y practicidad de la logística de riego, se sabe que dentro del lote hay un gradiente de pH en un intervalo de 3,5 a 5,5 pero no es posible realizar encalado suficiente para la corrección de este, ya que sale demasiado costoso y el agricultor no cuenta con los recursos económicos suficientes.
Factores:
R= Riego: Con 2 niveles (300 y 900 mm de lámina de riego). Este factor se asociará a las dos parcelas definidas por practicidad y factibilidad en campo, ya que es difícil controlar la logística y programación de riego para parcelas pequeñas, es mejor dejar que cada nivel de este factor esté en un área continua y lo más grande posible. Se definen estas dosis de riego con base a las necesidades del cultivo de forma que fueran contrastantes, ya que para un desarrollo óptimo de las plantas se reporta un requerimiento hídrico de 900 a 1000 mm por ciclo.
N= Fertilización Nitrogenada con urea: Con 4 niveles de dosis (50, 100,150 y 200 kg/ha). Este factor se asociará a las subparcelas, puesto que es más fácil manejar las dosis variables independientemente del tamaño del lote disponible para cada tratamiento.
Bloques:
De acuerdo con la información dada, se cuenta con una condición especial del suelo donde hay un gradiente de pH como se muestra en la gráfica siguiente, por lo cual se decide bloquear dejando en el Bloque 1 la zona del lote con pH más alto y en el Bloque 2 la zona con el pH más bajo.
Tratamientos: De acuerdo con las determinaciones anteriores los tratamientos quedan de la siguiente manera: N50R3 = Dosis de 50 kg/ha de N y 300 mm de riego N100R3 = Dosis de 100 kg/ha de N y 300 mm de riego N150R3 = Dosis de 150 kg/ha de N y 300 mm de riego N200R3 = Dosis de 200 kg/ha de N y 300 mm de riego N50R9 = Dosis de 50 kg/ha de N y 900 mm de riego N100R9 = Dosis de 100 kg/ha de N y 900 mm de riego N150R9 = Dosis de 150 kg/ha de N y 900 mm de riego N200R9 = Dosis de 200 kg/ha de N y 900 mm de riego
Cada tratamiento queda distribuido en ambos bloques para evaluar el efecto de estos.
Variable respuesta:
Se midió el rendimiento en g/planta de los tubérculos totales producidos y se muestrearon 10 plantas por tratamiento en la etapa de cosecha.
Covariable:
Dado que dentro del lote se presenta una variación importante de pH y gracias a que el agricultor cuenta con un medidor de pH portátil, se tiene pues la facilidad de realizar varias mediciones en diferentes partes del lote y de los bloques definidos para poder estimar si hay efecto de este parámetro sobre la variable evaluada o en interacción con los demás factores puesto que se parte del planteamiento que la disponibilidad de nutrientes está afectada por el pH, por lo cual se toma este parámetro como covariable ya que es un factor difícilmente controlable en el suelo.
Esquema del diseño Parcelas Divididas en Bloques Completos
Por medio del software R Studio, se realiza un diagrama de distribución y aleatorización del diseño planteado con ayuda de la función anm.ExpDesign para un diseño strip.split como el requerido para este experimento
anm.ExpDesign (method= "strip.split", titles= TRUE, cex.text= 1, mp.col= NULL, lwda= 1, n =8, EUcol= hcl.colors (n, palette="Dark 3"), interval= 0.5, iter= 2)
II) Analisis de covarianza
Revisión de supuestos
str(papa)## tibble [80 x 6] (S3: tbl_df/tbl/data.frame)
## $ Tratamiento : chr [1:80] "N50R3" "N50R3" "N50R3" "N50R3" ...
## $ Bloque : num [1:80] 1 1 1 1 1 1 1 1 1 1 ...
## $ Riego : num [1:80] 300 300 300 300 300 300 300 300 300 300 ...
## $ Fertilización: num [1:80] 50 50 50 50 50 50 50 50 50 50 ...
## $ Rendimiento : num [1:80] 345 491 373 118 275 ...
## $ pH : num [1:80] 4.42 4.43 4.52 4.55 4.56 ...papa$Riego =as.factor(papa$Riego)
papa$Fertilización =as.factor(papa$Fertilización)
str(papa)## tibble [80 x 6] (S3: tbl_df/tbl/data.frame)
## $ Tratamiento : chr [1:80] "N50R3" "N50R3" "N50R3" "N50R3" ...
## $ Bloque : num [1:80] 1 1 1 1 1 1 1 1 1 1 ...
## $ Riego : Factor w/ 2 levels "300","900": 1 1 1 1 1 1 1 1 1 1 ...
## $ Fertilización: Factor w/ 4 levels "50","100","150",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Rendimiento : num [1:80] 345 491 373 118 275 ...
## $ pH : num [1:80] 4.42 4.43 4.52 4.55 4.56 ...1) Homogeneidad de varianzas
leveneTest(papa$Rendimiento, papa$Fertilización, center="median") #Prueba Homocedasticidad factorA## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 3 1.7615 0.1617
## 76leveneTest(papa$Rendimiento, papa$Riego, center="median") #Prueba Homocedasticidad factorB## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 1 0.5404 0.4645
## 78Las varianzas son homogéneas
2) Independencia
anova1 = aov(papa$pH~papa$Fertilización) #anova covariable~factorA
summary(anova1)## Df Sum Sq Mean Sq F value Pr(>F)
## papa$Fertilización 3 0.667 0.2225 0.608 0.612
## Residuals 76 27.820 0.3660anova2 = aov(papa$pH~papa$Riego) #anova covariable~factorB
summary(anova2)## Df Sum Sq Mean Sq F value Pr(>F)
## papa$Riego 1 4.509 4.509 14.67 0.000258 ***
## Residuals 78 23.978 0.307
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Se cumple parcialmente el supuesto de independencia
3) Homogeneidad pendientes de regresión
anova3 = aov(papa$Rendimiento~papa$pH*papa$Riego)
summary(anova3)## Df Sum Sq Mean Sq F value Pr(>F)
## papa$pH 1 765348 765348 14.82 0.000245 ***
## papa$Riego 1 1072244 1072244 20.76 1.95e-05 ***
## papa$pH:papa$Riego 1 818541 818541 15.85 0.000156 ***
## Residuals 76 3925590 51653
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova4 = aov(papa$Rendimiento~papa$pH*papa$Fertilización)
summary(anova4)## Df Sum Sq Mean Sq F value Pr(>F)
## papa$pH 1 765348 765348 19.72 3.17e-05 ***
## papa$Fertilización 3 1681552 560517 14.45 1.84e-07 ***
## papa$pH:papa$Fertilización 3 1341048 447016 11.52 2.97e-06 ***
## Residuals 72 2793775 38802
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Revisión de supuestos para ANCOVA: Se encuentra que hay homogeneidad de varianzas e independencia del factor fertilización respecto a la covariable pH.
III) ANCOVA
ANCOVAt = aov(Rendimiento~pH+Tratamiento, data=papa) #efecto del tratamiento
summary(ANCOVAt)## Df Sum Sq Mean Sq F value Pr(>F)
## pH 1 765348 765348 19.59 3.40e-05 ***
## Tratamiento 7 3042555 434651 11.13 2.06e-09 ***
## Residuals 71 2773821 39068
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ANCOVAb = aov(Rendimiento~pH+Bloque, data=papa) #efecto del bloque
summary(ANCOVAb)## Df Sum Sq Mean Sq F value Pr(>F)
## pH 1 765348 765348 13.79 0.000385 ***
## Bloque 1 1542756 1542756 27.80 1.2e-06 ***
## Residuals 77 4273620 55502
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ANCOVAi = aov(Rendimiento~pH+Bloque:Tratamiento, data=papa) #interacción bloques:tratamientos
summary(ANCOVAi)## Df Sum Sq Mean Sq F value Pr(>F)
## pH 1 765348 765348 19.59 3.40e-05 ***
## Bloque:Tratamiento 7 3042555 434651 11.13 2.06e-09 ***
## Residuals 71 2773821 39068
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ANCOVAf = aov(Rendimiento~pH+Fertilización:Riego, data=papa) #interacción factores
summary(ANCOVAf)## Df Sum Sq Mean Sq F value Pr(>F)
## pH 1 765348 765348 19.59 3.40e-05 ***
## Fertilización:Riego 7 3042555 434651 11.13 2.06e-09 ***
## Residuals 71 2773821 39068
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1El análisis de covarianza muestra que los efectos del pH, de los tratamientos, del bloque y de la interacción bloque:tratamiento son altamente significativos para la variable rendimiento.
Análisis de Covarianza: Se encuentra que el pH, los bloques, los tratamientos, la interacción bloque:tratamiento y la interacción de los factores Fertilización:Riego tienen un efecto significativo sobre la variable respuesta, puesto que sus p valores son de 3.40e-05 para pH, 1.2e-06 para los bloques y 2.06e-09 para los tratamientos y las interacciones mencionadas. Entonces, se justifica el uso tanto de la covariable como del bloque en el modelo planteado
IV) Prueba de Comparación múltiple
TukeyHSD(ANCOVAt)## Warning in replications(paste("~", xx), data = mf): non-factors ignored: pH## Warning in TukeyHSD.aov(ANCOVAt): 'which' specified some non-factors which will
## be dropped## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Rendimiento ~ pH + Tratamiento, data = papa)
##
## $Tratamiento
## diff lwr upr p adj
## N100R9-N100R3 125.01728 -151.04820 401.082754 0.8476882
## N150R3-N100R3 291.08720 15.02172 567.152674 0.0315546
## N150R9-N100R3 335.69528 59.62980 611.760758 0.0070150
## N200R3-N100R3 -102.72442 -378.78990 173.341054 0.9398751
## N200R9-N100R3 355.39002 79.32454 631.455496 0.0034124
## N50R3-N100R3 -198.96545 -475.03093 77.100028 0.3351213
## N50R9-N100R3 23.01127 -253.05421 299.076744 0.9999955
## N150R3-N100R9 166.06992 -109.99556 442.135398 0.5697256
## N150R9-N100R9 210.67800 -65.38747 486.743482 0.2654976
## N200R3-N100R9 -227.74170 -503.80718 48.323778 0.1819625
## N200R9-N100R9 230.37274 -45.69273 506.438220 0.1709983
## N50R3-N100R9 -323.98273 -600.04820 -47.917248 0.0106024
## N50R9-N100R9 -102.00601 -378.07149 174.059469 0.9419910
## N150R9-N150R3 44.60808 -231.45739 320.673562 0.9996047
## N200R3-N150R3 -393.81162 -669.87710 -117.746142 0.0007696
## N200R9-N150R3 64.30282 -211.76266 340.368300 0.9958634
## N50R3-N150R3 -490.05265 -766.11812 -213.987168 0.0000128
## N50R9-N150R3 -268.07593 -544.14141 7.989548 0.0631964
## N200R3-N150R9 -438.41970 -714.48518 -162.354227 0.0001217
## N200R9-N150R9 19.69474 -256.37074 295.760216 0.9999985
## N50R3-N150R9 -534.66073 -810.72621 -258.595252 0.0000017
## N50R9-N150R9 -312.68401 -588.74949 -36.618536 0.0156065
## N200R9-N200R3 458.11444 182.04896 734.179920 0.0000522
## N50R3-N200R3 -96.24103 -372.30650 179.824452 0.9571577
## N50R9-N200R3 125.73569 -150.32979 401.801168 0.8438646
## N50R3-N200R9 -554.35547 -830.42095 -278.289991 0.0000007
## N50R9-N200R9 -332.37875 -608.44423 -56.313274 0.0078950
## N50R9-N50R3 221.97672 -54.08876 498.042194 0.2077625pairwise.t.test(papa$Rendimiento, papa$Tratamiento, p.adj = "bonf", paired=FALSE)##
## Pairwise comparisons using t tests with pooled SD
##
## data: papa$Rendimiento and papa$Tratamiento
##
## N100R3 N100R9 N150R3 N150R9 N200R3 N200R9 N50R3
## N100R9 0.00563 - - - - - -
## N150R3 0.11555 1.00000 - - - - -
## N150R9 1.1e-07 0.19403 0.01043 - - - -
## N200R3 1.00000 0.05470 0.75163 2.1e-06 - - -
## N200R9 0.00025 1.00000 1.00000 1.00000 0.00318 - -
## N50R3 1.00000 0.00026 0.00798 3.0e-09 1.00000 8.8e-06 -
## N50R9 1.00000 0.04216 0.61010 1.5e-06 1.00000 0.00236 1.00000
##
## P value adjustment method: bonferroniComparación de tratamientos: Prueba de comparación de medias usando Tukey.HSD y pairwise.t.test Con los resultados encontrados en la pairwise.t.test, donde se obtienen resultados más ajustados, se obtiene de manera general que hay diferencias significativas entre los tratamientos con riego de 900 mm respecto a los que tenían solo 300 mm y entre los tratamientos con mayor dosis de fertilización respecto a los que tenían la dosis más baja. Se resalta que el mejor tratamiento es N150R9, ya que es donde se obtiene el mayor rendimiento, presentando diferencias significativas con respecto a N200R3, esto se explica ya que el cultivo de papa es muy sensible al déficit hídrico y dado que tiene un requerimiento óptimo de 900 a 1000 mm por ciclo, claramente el tratamiento con menor dosis de riego afectó severamente el rendimiento, aunque tuviera la mayor dosis de fertilización con nitrógeno. Además también este mismo tratamiento presenta diferencias significativas respecto a N50R9, aunque este tenga buena disponibilidad de agua su dosis de fertilizante es insuficiente para lograr un adecuado crecimiento y desarrollo del cultivo, por tanto se comprueba la importancia que tiene la fertilización nitrogenada en la papa, sin embargo se tiene como recomendación hacer aplicación de N completamente en la siembra, ya que aplicar nitrógeno fraccionado puede ocasionar desbalance en la toma de potasio cuando la planta está en su estadio fenológico de tuberización y llenado, que es donde más necesita de este elemento (K). Se resalta también la importancia del pH en el suelo, ya que cuando se tienen valores muy bajos puede afectar la toma de nutrientes por parte de la planta, además de generar toxicidad si hay altos contenidos de aluminio en el suelo. Por otra parte, se plantea que la importancia en la interacción de los factores riego y fertilización, se debe a que para una adecuada toma de nutrientes es necesaria e indispensable la presencia de agua en el suelo, debido a que la mayor proporción de los nutrientes se mueven en flujos de masa con el agua hacia las raíces de la planta, lo que facilita su absorción y movimiento dentro de la planta.
Ejercicio Nº 8
I. Uso de los diseños en parcelas divididas
Cabe resaltar que el concepto de parcelas divididas se genera en el marco agrícola, por la necesidad de colocar en determinadas parcelas los factores que son difíciles de aleatorizar, (ejemplo: Riego, tipo de maquinaria, condiciones ambientales, variedades, fertilizante), esto no quiere decir que no se aleatoricen, sino que se hace de una manera estructurada y rígida, mientras que en las subparcelas los factores son fácilmente aleatorizados. En estos experimentos funciona como bloqueo uno de los factores, caso que no se presenta en los diseños completamente aleatorizados.
Cuando se plantean los diseños es muy importante identificar el nivel de aleatorización de cada factor, ya que suele confundirse en la práctica los diseños de parcelas divididas con los diseños completamente aleatorizados, permitiendo así que al momento de realizar el análisis de varianza no se contemplen los errores generados por la parcela, llegando a concluir erróneamente del resultado de los efectos obtenidos. Comparaciones entre estos dos diseños han demostrado que la eficiencia de dejar azar se reduce al 50%, lo que genera un sin sabor en cuanto al proceso de aleatorización.
Por lo general, los factores que se designan para las parcelas grandes tienden a ser muy costosos, por eso es recomendable realizar mas observaciones en las subparcelas y reducir en lo posible la aleatorización y la cantidad de parcelas completas, que ayudan a disminuir los costos y a su vez el tiempo del experimento.
Para determinar que diseño utilizar se deben tener en cuenta aspectos cualitativos y cuantitativos con el fin de cumplir con el objetivo del diseño, pero cabe destacar que un diseño de parcelas divididas bien planteado puede resultar practico y eficientes si se diseña correctamente, logrando así, que el investigador extraiga más información con un solo diseño.
Tomado de: [ASQ.org]
II. Sobre lo que significa unidad experimental y unidad de observación
A menudo suele confundirse estos dos términos debido a la concepción errónea que se tiene del concepto. La unidad experimental hace referencia a la unidad que recibe el tratamiento mientras que la unidad de observación es cada una de los individuos que son medidos u observados dentro de cada tratamiento, es decir, con un ejemplo muy simple, podríamos hablar de un experimento donde se tiene 1 parcela con papa, esta parcela a su vez se divide en 3 subparcelas, cada una tiene un tipo de fertilizante, como resultado tendríamos 3 tratamientos, de este diseño obtendríamos que la unidad experimental seria cada subparcela que tiene un fertilizante asignado, en total serian 3 unidades experimentales, pero resulta que dentro de cada subparcela se asignaron 30 plantas, entonces en este caso la unidad de observación seria cada planta que fue analizada, por lo tanto tendríamos un total de 90 unidades de observación. Es muy importante tener claro las diferencias entre estas dos unidades, ya que al incluir en el análisis de varianza el resultado obtenido por unidad de observación y no el experimental, tiene como consecuencia el aumento del error, lo que impacta al estadístico F obteniendo resultados erróneos.
Tomado de: [The Pennsylvania State University]
III. Guía para diseñar experimentos exitosos
Para iniciar es importante reconocer los cuatro pilares de un diseño experimental: Replicación, aleatorización, bloqueo y tamaño de las unidades experimentales. Por lo general los diseños requieren de presupuesto y tiempo, así que cada diseño debe estar enfocado “a que el fracaso no es una opción”, de lo cual se infiere que hay que reducir al máximo errores experimentales y tener un objetivo claro de lo que se quiere analizar, muchos de los fracasos se pueden evitar, por ejemplo, hay que evitar un mal diseño inicial, el mal planteamiento de los tratamientos y también la mala ejecución de los diseños y errores en la recopilación de los datos.
Los experimentos se dividen en tres tipos: Los observacionales, que analizan un comportamiento a través del tiempo sin que el investigador tenga algún control sobre los factores medidos o la población muestral. Los de medición, el investigador tiene control sobre el factor que se analiza, pero no de la muestra poblacional que a su vez no es aleatoria. Los comparables, donde el investigador tiene control de los factores y la población de muestreo. Para realizar diseños experimentales se tiene en cuenta los siguientes pasos:
- Se realizan preguntas e hipótesis respecto a un tema especifico
- Estas hipótesis se transforman en modelos del tema en cuestión
- Estos modelos se traducen en un modelo estadístico acompañado de un diseño estadístico, es aquí se identifican los tratamientos y el diseño a utilizar.
- Luego, se lleva a cabo el experimento con el cual se obtienen los resultados pertinentes
- Con los resultados se realiza un análisis estadístico y se interpretan los resultados obtenidos
- Por último, de acuerdo con los resultados obtenidos se generan nuevas hipótesis o modificaciones del diseño inicial.
Algo que rescata el autor son los dos tipos de información que se obtienen de los diseños experimentales, el análisis matemático o estadístico y las respuestas científicas enfocadas en el objeto de estudio, que en consecuencia brindan información para futuras investigaciones.
Replicación: Dentro de las funciones de este pilar se encuentra el estimar el error experimental a partir de varias repeticiones, también proporciona la precisión de un experimento, a su vez se aumenta el rango de condiciones encontradas durante el experimento y, por último, afecta el control del error. Aquí es muy importante definir las unidades experimentales y como replicar los tratamientos. Las replicaciones se pueden realizar a nivel de unidad experimental, replicación de todo el experimento de forma espacial y temporal, dentro de los niveles de la unidad experimental o medidas repetidas. En resumen, las réplicas tienen un impacto importante ya que permite detectar diferencias entre los tratamientos. En el campo de los diseños también es probable que no se realicen repeticiones, por lo tanto, los investigadores tienen como opción realizar experimentos utilizando bloques, también recurrir a diseños de parcelas de control, o una combinación de los dos anteriores. También es común utilizar diseños aumentados, donde se manejan un gran número de tratamiento, donde por lo general no se replica.
Aleatorización: Este pilar consiste en determinar si se realiza experimentos con factores parcialmente fijos (Parcelas divididas) o aleatorios (completamente al azar), con el fin de asignar al azar un tratamiento a la unidad experimental, este principio está ligado con el supuesto de independencia de los residuos. Este pilar tiene como función la estimación no sesgada de cada tratamiento y los errores experimentales y como medida de precaución contra cualquier perturbación que pueda surgir.
Bloqueo: Este pilar se utiliza con el fin de tener en cuenta los factores que puedan afectar la variable respuesta. Dentro de sus propósitos se encuentra crear grupos de unidades experimentales que sean homogéneos a diferencia de la aleatorización, o también para permitir diferentes tamaños de las unidades experimentales cuando se tienen parcelas mas grandes. Un error común es disponer de los bloques de forma lineal lo que genera errores en el análisis.
Tamaño de las unidades experimentales: Es uno de los principios menos estudiados pero que está ligado con la ley de Smith, donde se presenta una relación asintótica entre la varianza (por unida o parcela) y el tamaño de la parcela, donde a mayor tamaño de la parcela, menor es la varianza por unidad.
Como conclusiones, el no encontrar efectos significativos entre los tratamientos dan lugar a cuestionar el diseño, el procedimiento aplicado y los resultados obtenidos. Esto influye en la limitación para poder realizar publicaciones y acceder a bases de datos robustas acerca de investigaciones para poder mejorar la planificación de los diseños. Es de vital importancia tener claro los pilares del diseños de experimento y crear un diseño que responda tanto a los interrogantes y las hipótesis planteadas inicialmente como a los recursos y el objetivo que queremos estudiar, para esto tenemos una gran variedad de diseños que se adaptan a los diferentes tratamientos, bloques, aleatorizaciones, repeticiones, tamaño de muestras que se puedan presentar a la hora de planificar.
Tomado de: [Fundamentals of Experimental Design: Guidelines for Designing Successful Experiments]
Ejercicio Nº 9
Información del experimento
Variable respuesta: Rendimiento (kg/ha)
Factor I. Parcela principal, Arreglo espacial: Distancia entre surcos (0.7 y 0.8 metros)
Factor II. Subparcelas, Arreglo espacial: Numero de plantas por metro lineal (7,8 y 9 plantas)
Tratamientos: 6
Repeticiones: 4
a. Estructura factorial: Diseño en bloques completos al azar en parcelas divididas.
b. La razón de colocar cada factor en la parcela principal o en la subparcela: En la parcela principal se fijó la distancia entre hileras, esto se debe a que, comparando con el segundo factor, es más practico por temas de sembrado, ya sea utilizando maquinas o de forma manual, definir el mismo ancho de surco y no estar variando entre cada hilera. Mientras que el número de plantas, que para este caso se realizó de manera manual, resultaba más fácil aleatorizarlo en las subparcelas.
c. La revisión de supuestos para el análisis de varianza: De acuerdo con lo que reporta el autor, se realizó un análisis preliminar en cada uno de los análisis univariado o multivariado, pero no se reporta resultado alguno sobre los supuestos en el artículo. En este caso nos veremos obligados a creer que lo que dice el autor es cierto y que para poder concluir respecto a cada variable se cumplieron con dichos supuestos.
d. La tabla del análisis de varianza: En el artículo como tal, no se muestra el análisis de varianza realizado por cada variable con el software o algo similar, sin embargo, se adjuntaron algunas tablas-resumen donde se especificaba donde habían diferencias significativas entre los factores, la interacción de los factores y cada una de las variables medidas durante el experimento cuando se realizó el análisis univariado; mientras que para los análisis multivariados no se reportaron tablas con este resumen, solo suministraron tablas comparando medias de cada variable evaluada. El autor recalca en cada gráfica y tabla que la asignación de “letras” determina si hubo o no diferencias significativas, haciendo referencia a las pruebas de comparación de medias.
e. El uso de muchos análisis de varianzas en lugar de uno solo multivariante: Este investigador uso los dos tipos de análisis, por ejemplo, cuando observaba el efecto de los factores y la interacción de los factores en variables como: Clorofila, altura de la planta, altura de la mazorca y área del tallo, utilizo el análisis de varianza univariado; mientras para otro tipo de variables como: Características morfológicas (Índice de semillas, Diámetro de la mazorca, Tamaño de la mazorca, Numero de filas en la mazorca, numero de granos por fila y numero de granos por mazorca), Índice de área foliar, coeficiente de extinción de luz y radiación solar interceptada por día, se empleó el análisis multivariado de varianza. Para el caso del rendimiento, no reporta el tipo de análisis realizado. Desde nuestro punto de vista, el análisis parece ser el correcto dada la complejidad de evaluar la correlación entre las diferentes variables independientes y dependientes, para poder concluir con un análisis más completo.
f. El método de comparaciones de medias después del Anova: Para el caso de análisis univariado se realizó la comparación de medias con la prueba de Tukey, mientras que para el análisis multivariado se utilizó la prueba de Hotelling.
g. La interacción de factores: De acuerdo con el autor, en el análisis de los resultados indica que la única interacción entre la distancia de surcos y el número de plantas fue significativa cuando se analizó el rendimiento, tanto por kg/ha como de g/planta. Para el primer caso, se cumple que, a menor distancia entre hilera y mayor número de plantas por metro líneal, aumenta el rendimiento, ya que hay mayor densidad de siembra. Mientras que, en el segundo caso, el mejor rendimiento por planta se encontró en la mayor distancia entre surco y menor número de plantas, dado que las plantas no tenían que competir por los recursos presentes en el suelo, recibían mayor extensión de luz y las condiciones eran menos estresantes para la planta.
h. El balanceo o desbalanceo: No se menciona nada al respecto en el artículo, son explicitos en decir que se realizaron las repeticiones correspondientes pero de ahí a decir que hizo falta algun dato creo que no sería benefico en lo absoluto para el autor del articulo.
i. La definición clara de la unidad experimental: El autor es claro en decir como son los tratamientos y al señalar que en cada parcela (unidad experimental) donde se utiliza el factor de distancia entre surcos se colocaron subparcelas aleatorizadas para el número de plantas por hilera por cada metro lineal.
j. Software utilizado y librería específica: El software utilizado fue R, en cuanto a las libreriras no se hace mención alguna en el articulo.