Untitled
Sara Daniela Vanegas Garzón
10/11/2020
Código:
U: 3
T: 1
Nota: se trabajan estos ejercicios con un nivel de confianza del 95%.
EJERCICIO 1
a-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 8U3. Con esta información complete la tabla del ANOVA que se muestra a continuación. Si el F tabulado es 2.8. ¿qué puede decirse acerca de la Hipótesis nula de igualdad de los promedios del índice en todas las condiciones de tratamiento (use el p valor sí como el cociente F calculado de la tabla para concluir)? ¿vale la pena comparar las medias de tratamientos a posteriori del ANOVA (prueba de Tukey)? b-De ser afirmativa su respuesta aplique este procedimiento y concluya (seleccione el o los tratamientos con mayores contenidos de clorofila). Para esto genere unos datos para cada tratamiento que tengan una varianza como la dada en el ejercicio y una media que desde el primer al último tratamiento cambie a razón de 10 unidades del índice, empezando con 40U. data1=rnorm(8,405,sqrt(853)) ;data2=rnorm(8,415,sqrt(853));…;data9= Revise los tres supuestos del ANOVA discutidos en clase y con este valide el uso del ANOVA en este estudio (con los datos generados). De no cumplirse el supuesto de homocedasticidad utilice el ANOVA de Welch discutido en clase y vuelva a concluir.
a. Tabla del análisis de varianza
Datos que se tienen:
Número de datos obtenidos: 72
Número de tratamientos: 9
Varianza de las 72 observaciones: 833 = SC tot/GL tot
Número de repeticiones: 8
Se sabe que:
La suma de cuadrados (SCt) = (SC trt + SC rep)
Grados de libertad tratamientos: 9-1 = 8
Grados de libertad totales: 72-1 =71
Grados de libertad repeticiones: (71-8) = 63
Varianza de los tratamientos: SC trt/Gl trt
Varianza de las repeticiones: SC rep/Gl rep
Resolviendo las ecuaciones con los datos se obtiene la siguiente tabla:
Tabla del análisis de varianza (ANOVA)
El \(F\) calculado dio \(0.8891\), al ser tan pequeño quiere decir que las diferencias son causadas mayormente por las repeticiones que por los tratamientos, esto no es deseable en los experimentos agronómicos principalmente en este caso.
Planteo de hipótesis para los tratameintos
\[Ho: \mu_1= \mu_2= \mu_3 =\mu_4 =\mu_5= \mu_6= \mu_7 =\mu_8=\mu_9 \\Ha: Ho \ es \ falsa\]
PrF=pf(q=0.89, df1=8, df2=63, lower.tail = F); PrF #P-valor## [1] 0.5301822Curva F
d=df(seq(0,5,0.1),8,63)
plot(seq(0,5,0.1), d, type="l", main = "Curva F")
ftab= 2.8; ftab## [1] 2.8abline(v=ftab, col="red")
text(1,0.2, "No rechazo Ho")
text(4,0.2, "Rechazo Ho")
text(3,0.02, "5%")
arrows(x0 = 0.8891, y0= 0.4, x1=0.8891, y1=0, col="blue")
text(0.8891,0.45, "Fc=0.8891")
text(ftab, 0.5, "Ftab")
Desde un punto de vista estadístico si el \(F\) tabulado es de 2,8 indica que NO se rechaza la hipótesis nula de igualdad de promedios de los tratamientos debido a que el \(F\) calculado (0.8891) dio menor al \(F\) tabulado. Por tanto, los tratamientos no difieren estadísticamente hablando.
- No vale la pena comparar las medias de los tratamientos con la prueba de Tukey, a pesar de que esta compara las medias de los tratamientos, debido a que en este ejemplo las repeticiones fueron los que causaron las diferencias más no los tratamientos, por tanto no es conveniente utilizarla.Por otro lado, como no se rechazó la hipótesis nula eso quiere decir que la media de los tratamientos se comportan IGUAL de ahí que, se escoja cual se escoja van a dar resultados similares, esto estadísticamente hablando.
EJERCICIO 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 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.
a-Realice el ANOVA para este diseño recordando que es un caso desbalanceado. Concluya sobre el resultado de la tabla del ANOVA obtenida. (¿Afecta el orden de colocación de los efectos del modelo dentro del software R? Verifique si la tabla del ANOVA cambia): en R: lm() lm(respuesta~ bloque+método) lm(respuesta~ método + bloque) b-Estimar el valor de la observación usando el promedio de los datos para los cinco granjeros del mismo método M1 y luego realice el análisis de varianza para probar las diferencias en las pérdidas medias de peso para los cinco métodos de limpiado de las fibras de algodón. Compare este resultado con el caso desbalanceado (de ser posible).
Modelo del diseño con sus condiciones laterales
\[Y_{ij} = \mu + \tau_i + \beta_j + \epsilon_{ij}\\ i= 1\dots5;\ j=1\dots6 \\ \sum_{i=1}^5 \tau_i=0 \\ \sum_{i=1}^6 \beta_j=0\]
Modelo de 2 vías, se tiene el bloque granja y se tiene el factor método, pero solo hay 1 observación por bloque.
# Este diseño es factorial simple en bloques desbalanceado debido a que hace falta un dato, usé NA que reemplaza el dato faltante.
perdidas = c(NA,6.75,13.05,10.26,8.01,8.42,5.54,3.53,11.20,7.21,3.24,6.45,7.67,4.15,9.79,8.27,6.75,5.50,7.89,1.97,8.97,6.12,4.22,7.84,9.27,4.39,13.44,9.13,9.20,7.13)
perdidas## [1] NA 6.75 13.05 10.26 8.01 8.42 5.54 3.53 11.20 7.21 3.24 6.45
## [13] 7.67 4.15 9.79 8.27 6.75 5.50 7.89 1.97 8.97 6.12 4.22 7.84
## [25] 9.27 4.39 13.44 9.13 9.20 7.13granja = c("G1","G2","G3","G4","G5","G6","G1","G2","G3","G4","G5","G6","G1","G2","G3","G4","G5","G6","G1","G2","G3","G4","G5","G6","G1","G2","G3","G4","G5","G6"); granja #actua como bloque, es una razon de bloqueo## [1] "G1" "G2" "G3" "G4" "G5" "G6" "G1" "G2" "G3" "G4" "G5" "G6" "G1" "G2" "G3"
## [16] "G4" "G5" "G6" "G1" "G2" "G3" "G4" "G5" "G6" "G1" "G2" "G3" "G4" "G5" "G6"metodo = gl(5,6,30, labels=c("M1","M2","M3","M4","M5"))
metodo## [1] M1 M1 M1 M1 M1 M1 M2 M2 M2 M2 M2 M2 M3 M3 M3 M3 M3 M3 M4 M4 M4 M4 M4 M4 M5
## [26] M5 M5 M5 M5 M5
## Levels: M1 M2 M3 M4 M5df2 = data.frame(perdidas, granja, metodo)
df2## perdidas granja metodo
## 1 NA G1 M1
## 2 6.75 G2 M1
## 3 13.05 G3 M1
## 4 10.26 G4 M1
## 5 8.01 G5 M1
## 6 8.42 G6 M1
## 7 5.54 G1 M2
## 8 3.53 G2 M2
## 9 11.20 G3 M2
## 10 7.21 G4 M2
## 11 3.24 G5 M2
## 12 6.45 G6 M2
## 13 7.67 G1 M3
## 14 4.15 G2 M3
## 15 9.79 G3 M3
## 16 8.27 G4 M3
## 17 6.75 G5 M3
## 18 5.50 G6 M3
## 19 7.89 G1 M4
## 20 1.97 G2 M4
## 21 8.97 G3 M4
## 22 6.12 G4 M4
## 23 4.22 G5 M4
## 24 7.84 G6 M4
## 25 9.27 G1 M5
## 26 4.39 G2 M5
## 27 13.44 G3 M5
## 28 9.13 G4 M5
## 29 9.20 G5 M5
## 30 7.13 G6 M5M1=c(NA,6.75,13.05,10.26,8.01,8.42) # datos del metodo 1
M1## [1] NA 6.75 13.05 10.26 8.01 8.42mediaM1=mean(na.omit(M1)) # media del M1 sin contar el dato faltante
mediaM1## [1] 9.298a.
HIPÓTESIS PARA EL FACTOR
\[Ho: \mu_1= \mu_2= \mu_3 =\mu_4 =\mu_5 \\Ha: Ho \ es \ falsa\]
library(car) ## Loading required package: carData# usando la libreria car
# Primero bloque y luego factor, este es el método correcto, con anova
mod1b= lm(perdidas~granja+metodo, data = df2) # uso de la funcion lm()
anova(mod1b)## Analysis of Variance Table
##
## Response: perdidas
## Df Sum Sq Mean Sq F value Pr(>F)
## granja 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 ' ' 1library(car)
# usando la libreria car
# Primero bloque y luego factor, este es el método correcto, con Anova
mod1c= lm(perdidas~granja+metodo, data = df2) # uso de la funcion lm()
Anova(mod1c)## Anova Table (Type II tests)
##
## Response: perdidas
## Sum Sq Df F value Pr(>F)
## granja 139.66 5 20.2331 5.163e-07 ***
## metodo 49.12 4 8.8951 0.0003186 ***
## Residuals 26.23 19
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1library(car) # usando la libreria car
# Ahora primero factor y luego bloque, esta es forma no es correcta
mod1d= Anova(lm(perdidas~metodo+granja, data = df2)) # uso de la funcion lm()
mod1d## Anova Table (Type II tests)
##
## Response: perdidas
## Sum Sq Df F value Pr(>F)
## metodo 49.12 4 8.8951 0.0003186 ***
## granja 139.66 5 20.2331 5.163e-07 ***
## Residuals 26.23 19
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1library(car) # usando la libreria car
# Ahora primero factor y luego bloque, esta es forma no es correcta
mod1e= lm(perdidas~metodo+granja, data = df2) # uso de la funcion lm()
anova(mod1e)## Analysis of Variance Table
##
## Response: perdidas
## Df Sum Sq Mean Sq F value Pr(>F)
## metodo 4 47.763 11.9407 8.6494 0.0003754 ***
## granja 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 ' ' 1Se observa que los valores son similares entre los dos análisis corridos todos indican p.valores significativos Para correr este diseño desbalanceado se usó la función lm(), la cual sirve para correr diseños desbalanceados. Como se puede observar con los análisis de varianza presentados dan prácticamente igual si se corre primero con el factor método y luego bloque granja que si se corre primero el bloque y luego el factor, eso sí, usando la función lm(). Sin embargo, lo más recomendable es colocar primero el bloque y luego si el factor que se está analizando.
Los resultados del ANOVA indican es que el p-valor del método (0.0003186) dio significativo, menor al nivel de significancia del 5% o 0.05, lo cual se concluye que existen diferencias estadísticas en los métodos de limpieza de las fibras de algodón. Se rechaza la hipótesis nula de que las medias de los 5 tratamientos o métodos son iguales, esto significa que al menos una media es diferente. Esto teniendo en cuenta que hizo falta un dato en el método 1.
library(lattice) # usando la libreria lattice
dotplot(perdidas~metodo|granja) #grafico correcto porque los metodos estan dentro de las granjas, esta es la realidad
#dotplot(perdidas~granja|metodo) #grafico incorrecto Los tratamientos o métodos en este caso están dentro de cada granja
b.
Imputación con la media para que ahora el diseño sea balanceado, acá se tomó el promedio de los 5 granjeros del mismo método 1
perdidasM = c(9.298 , 6.75,13.05,10.26,8.01,8.42,5.54,3.53,11.20,7.21,3.24,6.45,7.67,4.15,9.79,8.27,6.75,5.50,7.89,1.97,8.97,6.12,4.22,7.84,9.27,4.39,13.44,9.13,9.20,7.13)
perdidasM## [1] 9.298 6.750 13.050 10.260 8.010 8.420 5.540 3.530 11.200 7.210
## [11] 3.240 6.450 7.670 4.150 9.790 8.270 6.750 5.500 7.890 1.970
## [21] 8.970 6.120 4.220 7.840 9.270 4.390 13.440 9.130 9.200 7.130granja = c("G1","G2","G3","G4","G5","G6","G1","G2","G3","G4","G5","G6","G1","G2","G3","G4","G5","G6","G1","G2","G3","G4","G5","G6","G1","G2","G3","G4","G5","G6"); granja #actua como bloque## [1] "G1" "G2" "G3" "G4" "G5" "G6" "G1" "G2" "G3" "G4" "G5" "G6" "G1" "G2" "G3"
## [16] "G4" "G5" "G6" "G1" "G2" "G3" "G4" "G5" "G6" "G1" "G2" "G3" "G4" "G5" "G6"metodo = gl(5,6,30, labels=c("M1","M2","M3","M4","M5"))
metodo## [1] M1 M1 M1 M1 M1 M1 M2 M2 M2 M2 M2 M2 M3 M3 M3 M3 M3 M3 M4 M4 M4 M4 M4 M4 M5
## [26] M5 M5 M5 M5 M5
## Levels: M1 M2 M3 M4 M5df3 = data.frame(perdidasM, granja, metodo)
df3## perdidasM granja metodo
## 1 9.298 G1 M1
## 2 6.750 G2 M1
## 3 13.050 G3 M1
## 4 10.260 G4 M1
## 5 8.010 G5 M1
## 6 8.420 G6 M1
## 7 5.540 G1 M2
## 8 3.530 G2 M2
## 9 11.200 G3 M2
## 10 7.210 G4 M2
## 11 3.240 G5 M2
## 12 6.450 G6 M2
## 13 7.670 G1 M3
## 14 4.150 G2 M3
## 15 9.790 G3 M3
## 16 8.270 G4 M3
## 17 6.750 G5 M3
## 18 5.500 G6 M3
## 19 7.890 G1 M4
## 20 1.970 G2 M4
## 21 8.970 G3 M4
## 22 6.120 G4 M4
## 23 4.220 G5 M4
## 24 7.840 G6 M4
## 25 9.270 G1 M5
## 26 4.390 G2 M5
## 27 13.440 G3 M5
## 28 9.130 G4 M5
## 29 9.200 G5 M5
## 30 7.130 G6 M5# Ahora corriendo el análisis balanceado
mod3= aov(perdidasM~granja+metodo, data=df3)
summary(mod3)## Df Sum Sq Mean Sq F value Pr(>F)
## granja 5 139.36 27.873 21.015 2.45e-07 ***
## metodo 4 51.15 12.787 9.641 0.000164 ***
## Residuals 20 26.53 1.326
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Como se observa en el ANOVA existen diferencias estadísticas en las pérdidas medias de peso de los cinco métodos de limpiado de las fibras del algodón, debido a que su p.valor fue significativo (0.000164) dio menor al 0.05. Por tanto se rechaza la hipóteisis nula, de ahí que al menos uno es diferente. El p.valor del bloque no se le hace un análisis, sin embargo este al dar tan pequeño muestra que si valió la pena bloquear por granja. Al comparar este ANOVA con el ANOVA anterior los dos llegan a la misma respuesta: se rechaza la hipótesis nula y por tanto los tratamientos son diferentes.
tk=TukeyHSD(mod3)
tk## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = perdidasM ~ granja + metodo, data = df3)
##
## $granja
## diff lwr upr p adj
## G2-G1 -3.7756 -6.0650872 -1.4861128 0.0005620
## G3-G1 3.3564 1.0669128 5.6458872 0.0020390
## G4-G1 0.2644 -2.0250872 2.5538872 0.9990388
## G5-G1 -1.6496 -3.9390872 0.6398872 0.2535915
## G6-G1 -0.8656 -3.1550872 1.4238872 0.8371590
## G3-G2 7.1320 4.8425128 9.4214872 0.0000001
## G4-G2 4.0400 1.7505128 6.3294872 0.0002515
## G5-G2 2.1260 -0.1634872 4.4154872 0.0784089
## G6-G2 2.9100 0.6205128 5.1994872 0.0080275
## G4-G3 -3.0920 -5.3814872 -0.8025128 0.0046011
## G5-G3 -5.0060 -7.2954872 -2.7165128 0.0000148
## G6-G3 -4.2220 -6.5114872 -1.9325128 0.0001454
## G5-G4 -1.9140 -4.2034872 0.3754872 0.1361201
## G6-G4 -1.1300 -3.4194872 1.1594872 0.6372884
## G6-G5 0.7840 -1.5054872 3.0734872 0.8850199
##
## $metodo
## diff lwr upr p adj
## M2-M1 -3.10300000 -5.0926820 -1.1133180 0.0012505
## M3-M1 -2.27633333 -4.2660153 -0.2866514 0.0201403
## M4-M1 -3.12966667 -5.1193486 -1.1399847 0.0011421
## M5-M1 -0.53800000 -2.5276820 1.4516820 0.9247537
## M3-M2 0.82666667 -1.1630153 2.8163486 0.7270386
## M4-M2 -0.02666667 -2.0163486 1.9630153 0.9999994
## M5-M2 2.56500000 0.5753180 4.5546820 0.0077482
## M4-M3 -0.85333333 -2.8430153 1.1363486 0.7039395
## M5-M3 1.73833333 -0.2513486 3.7280153 0.1053894
## M5-M4 2.59166667 0.6019847 4.5813486 0.0070847# medias de las 5 granjas
mediasM=tapply(df3$perdidasM,c(df3$granja), mean)
mediasM## G1 G2 G3 G4 G5 G6
## 7.9336 4.1580 11.2900 8.1980 6.2840 7.0680# medias de los 5 metodos
mediasM1=tapply(df3$perdidasM,c(df3$metodo), mean)
mediasM1## 1 2 3 4 5
## 9.298000 6.195000 7.021667 6.168333 8.760000library(lattice)
dotplot(perdidasM~metodo|granja)
Al final se observa que la granja 2 tuvo menores pérdidas en los 5 métodos, en esta el método 4 fue el que presentó el menor dato de pérdida de algodón (1.97 kg). Sin embargo, lo que interesa en este estudio es saber cual método conviene más para evitar menos daños y pérdidas en las fibras del algodón, por tanto se recomienda el uso del método 4 y 2 que tuvieron menos pérdidas, un promedio de 6.17 y 6.20 respectivamente, ya va depender del investigador de escoger el método que prefiera ya sea mecánico o una combinación mecánica y química.
EJERCICIO 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.U+0.1 y de la capa inferior osciló entre 2 y 2.T+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+U 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%.
Generación de datos
set.seed(2021) # se fija semilla
CO_5 = sort.int(runif(50,3.0, 3.4), partial = 28) #5 cm de profundidad
CO_5## [1] 3.089954 3.096026 3.069193 3.152698 3.013110 3.170538 3.054168 3.106672
## [9] 3.107394 3.156998 3.010907 3.006674 3.055832 3.020835 3.100629 3.180507
## [17] 3.202198 3.220906 3.226981 3.218279 3.241296 3.178934 3.254530 3.280538
## [25] 3.228883 3.256176 3.283873 3.313512 3.317543 3.325996 3.329332 3.338663
## [33] 3.334996 3.326169 3.347015 3.328021 3.372285 3.360837 3.373960 3.379095
## [41] 3.377084 3.378278 3.352194 3.382922 3.373717 3.364858 3.354599 3.382136
## [49] 3.383273 3.393195CO_10 = sort.int(runif(50,2,2.5), partial = 28) #10 cm de profundidad
CO_10## [1] 2.117028 2.107068 2.036387 2.108672 2.144516 2.093710 2.036439 2.049617
## [9] 2.028328 2.023944 2.012388 2.136372 2.036908 2.038799 2.058136 2.180832
## [17] 2.198033 2.182502 2.198839 2.190546 2.243144 2.261838 2.204301 2.261845
## [25] 2.231337 2.227485 2.232654 2.265601 2.270981 2.311312 2.273164 2.340401
## [33] 2.346564 2.331138 2.339410 2.328727 2.334111 2.347641 2.410992 2.364244
## [41] 2.377760 2.381717 2.464981 2.491498 2.480407 2.448985 2.497622 2.423002
## [49] 2.468699 2.432240#Creacion de la ventana
VEN=expand.grid(longitud=seq(0,100,25), latitud=seq(0,200,length.out = 10))
head(VEN)## longitud latitud
## 1 0 0.00000
## 2 25 0.00000
## 3 50 0.00000
## 4 75 0.00000
## 5 100 0.00000
## 6 0 22.22222gg=data.frame(Longitud=rep(VEN$longitud,2),Latitud=rep(VEN$latitud,2), Profundidad = rep(c(-5,-10), each = 50), col =c(CO_5,CO_10))
head(gg)## Longitud Latitud Profundidad col
## 1 0 0.00000 -5 3.089954
## 2 25 0.00000 -5 3.096026
## 3 50 0.00000 -5 3.069193
## 4 75 0.00000 -5 3.152698
## 5 100 0.00000 -5 3.013110
## 6 0 22.22222 -5 3.170538Punto de vista 3D
library(plotly)## Loading required package: ggplot2##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutplot_ly(x=gg$Latitud, y=gg$Longitud, z=gg$Profundidad, type="scatter3d", mode="markers", color = gg$col)## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.Se puede observar a simple vista que hay mas contenido de carbono organico a una profundidad de 5 cm (-5 cm) que a una profundidad de 10 cm (-10 cm). Por otro lado, tiende a variar mas el carbono organico en la latitud, a medida que se avanza en latitud (de 0 a 200) el C.O tiende a aumentar, esto se puede observar en las dos capas. Mientras que en la longitud parece constante mas o menos el valor de C.O. Los mayores contenidos de C.O se encuentran en la latitud 200, longotud 100 a una profundidad del suelo de 5 cm.
# explorando la correlacion
#cor(CO_5, CO_10)
#plot(CO_5, CO_10, xlim =c(3, 3.3), ylim =c(2,2.3), pch=16, main = "Correlacion") # Parece existir correlacion entre los datosComo los datos provienen de una distribución uniforme, es probable que sus promedios no sean normales, por lo tanto se va verificar si son o no normales para asi luego decidir que prueba aplicar (si t.test o wilcox.test) al momento de comparar las capas.
Comparación de capas
Verificando la normalidad de los promedios de cada capa
\[H_o: Promedios \ normales \\ H_a: Promedios \ NO \ normales\]
# Uso shapiro.test que indica si son o no normales los promedios de la capa superior
test_norm= shapiro.test(CO_5); test_norm##
## Shapiro-Wilk normality test
##
## data: CO_5
## W = 0.89479, p-value = 0.0003247Se puede concluir que como el p.valor dio menor a 0.05, se rechaza la hipotesis nula y por tanto los promedios de la capa superior no son normales.
Prueba de Shapiro busqueda de normalidad
# Uso shapiro.test que indica si son o no normales los promedios de la capa inferior
test_norm1= shapiro.test(CO_10); test_norm1##
## Shapiro-Wilk normality test
##
## data: CO_10
## W = 0.95175, p-value = 0.04028Con un nivel de confianza del 95% se puede concluir que como el p.valor dio menor a 0.05, se rechaza la hipotesis nula y por tanto los promedios de la capa inferior tampoco son normales. Se concluyó que los promedios no siguen una distribución normal.
Prueba de varianzas
\[H_o: \sigma^2_{CO_5} = \sigma^2_{CO_{10}} \\ H_a: \sigma^2_{CO_5} \neq \sigma^2_{CO_{10}}\]
prueba_var=var.test(CO_5,CO_10); prueba_var##
## F test to compare two variances
##
## data: CO_5 and CO_10
## F = 0.7254, num df = 49, denom df = 49, p-value = 0.2646
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.4116477 1.2782930
## sample estimates:
## ratio of variances
## 0.7254008Se concluye que las varianzas son iguales
Ahora se tendrá que utilizar la prueba de Wilcoxon para comparar las capas
\[H_o: \ Med_{CO_5} \ = \ Med_{CO_{10}} \\ H_a: \ Med_{CO_5} \neq \ Med_{CO_{10}}\]
Prueba de Wilcoxon para dos muestras pareadas espacialmente
prbw = wilcox.test(CO_5, CO_10, paired = TRUE, mu = 0, conf.level = 0.95, alternative = 't'); prbw##
## Wilcoxon signed rank test with continuity correction
##
## data: CO_5 and CO_10
## V = 1275, p-value = 7.79e-10
## alternative hypothesis: true location shift is not equal to 0ifelse(prbw$p.value<0.05,'Rechazo Ho','No rechazo Ho')## [1] "Rechazo Ho"Como se pudo observar y utilizando un nivel de confianza del 95% se rechaza la hipotesis nula, ya que el p.valor (7.79e-10) dio menor a 0.05. se rechaza la hipotesis nula de que la media de carbono organico a 5 y 10 cm son iguales. Por tanto si se encuentran diferencias en la media de carbono organico entre las capas a profundidades de 5 y 10 cm.
EJERCICIO 4.
El siguiente diseño se corresponde con un factorial completo (32) en arreglo completamente al azar. Los factores y la respuesta fueron creados con el código:
• Escriba (completamente especificado) el modelo del diseño • Realice el Anova para este diseño y de ser necesario realice la prueba de comparaciones de medias para los efectos principales de F1: dosis de un insecticida que se sospecha tiene un efecto de disminución del crecimiento (biomasa) y F2: número de aplicaciones durante el desarrollo del cultivo. • Use los resultados del ANOVA y el gráfico de interacción (ggplot2) para visualizar si existe o no interacción entre los factores. El investigador quiso colocar como covariable el contenido de arcilla(expansible) en el suelo utilizado en cada unidad experimental. Genere unos datos con la distribución uniforme cuya medida oscile entre 0.20 y 0.40 , ordene estas medidas en forma decreciente y meta dentro del análisis esta variable. Especifique nuevamente el modelo y realice el análisis de covarianza respectivo ¿se justifica el uso de la covariable? Construya nuevamente el gráfico de interacción y compare con el caso sin covariable (discuta el resultado). Revise en internet los supuestos que deben tener las covariables para ser utilizadas en el modelo. ¿Se está incumpliendo en nuestros datos alguno de los supuestos necesarios? Revise los supuestos sobre los residuales tanto del ANOVA como del ANCOVA ¿qué puede percibir? ¿recomendaría el uso de arcillas para minimizar el efecto sobre el contenido de biomasa que puede ocasionar el uso del insecticida?
set.seed(2020) # fijar semilla
D=expand.grid(F1 = c(3.25,3.75,4.25), F2 = c(4,5,6)); D # Crea el diseNo 3^2## F1 F2
## 1 3.25 4
## 2 3.75 4
## 3 4.25 4
## 4 3.25 5
## 5 3.75 5
## 6 4.25 5
## 7 3.25 6
## 8 3.75 6
## 9 4.25 6D=rbind(D,D); D # Crea la estructura para dos repeticiones por tratamiento## F1 F2
## 1 3.25 4
## 2 3.75 4
## 3 4.25 4
## 4 3.25 5
## 5 3.75 5
## 6 4.25 5
## 7 3.25 6
## 8 3.75 6
## 9 4.25 6
## 10 3.25 4
## 11 3.75 4
## 12 4.25 4
## 13 3.25 5
## 14 3.75 5
## 15 4.25 5
## 16 3.25 6
## 17 3.75 6
## 18 4.25 6D=D[order(sample(1:18)),]; D # aleatoriza la estructura## F1 F2
## 2 3.75 4
## 10 3.25 4
## 16 3.25 6
## 4 3.25 5
## 13 3.25 5
## 6 4.25 5
## 17 3.75 6
## 8 3.75 6
## 14 3.75 5
## 5 3.75 5
## 9 4.25 6
## 1 3.25 4
## 12 4.25 4
## 11 3.75 4
## 15 4.25 5
## 7 3.25 6
## 3 4.25 4
## 18 4.25 6D$biomasa=sort.int(rnorm(18,3,0.3), partial=9); D$biomasa# Crea la respuesta## [1] 2.708826 2.772692 2.143359 2.560519 2.708666 2.773705 2.770350 2.832470
## [9] 2.898280 3.359619 3.054099 3.157896 3.487669 3.451547 3.016111 3.042156
## [17] 3.200552 2.989329#F1 es el factor dosis de un insecticida
#F2 es el factor numero de aplicacionesModelo del diseño factorial completo completamente al azar
\[y_{ijk} = \mu + \tau_i + \beta_j + (\tau\beta)_{ij} + \epsilon_{ijk}\\ i=1,2,3,\\ j=1,2,3;\\ k=1,2 \\ \ y \ sus \ condiciones \ laterales \\ \sum_{i=1}^5 \tau_i=0 \\ \sum_{i=1}^6 \beta_j=0\]
\[H_{01}:(\tau\beta)_{ij}=0;\forall_{ij}\]
\[H_{02}:\tau_i=0;\forall_{i}\] \[H_{03}:\beta_j=0;\forall_{j}\]
ANOVA
# ANOVA
D$F1=as.factor(D$F1); D$F1## [1] 3.75 3.25 3.25 3.25 3.25 4.25 3.75 3.75 3.75 3.75 4.25 3.25 4.25 3.75 4.25
## [16] 3.25 4.25 4.25
## Levels: 3.25 3.75 4.25D$F2=as.factor(D$F2); D$F2## [1] 4 4 6 5 5 5 6 6 5 5 6 4 4 4 5 6 4 6
## Levels: 4 5 6mod=aov(D$biomasa ~ F1 * F2, data= D)
summary(mod)## Df Sum Sq Mean Sq F value Pr(>F)
## F1 2 0.4161 0.20804 1.979 0.194
## F2 2 0.3426 0.17128 1.630 0.249
## F1:F2 4 0.1635 0.04087 0.389 0.812
## Residuals 9 0.9459 0.10510Primero: el p.valor de la interaccion dio mayor al nivel de significancia del 5% o 0.05, lo que significa que NO se rechaza la hipotesis nula indicando que el modelo no muestra interacciOn (ver mas adelante lo que se concluyo).
# Corriendo el modelo sin interaccion
mod9=aov(D$biomasa ~ F1 + F2, data=D)
summary(mod9)## Df Sum Sq Mean Sq F value Pr(>F)
## F1 2 0.4161 0.20804 2.438 0.126
## F2 2 0.3426 0.17128 2.007 0.174
## Residuals 13 1.1094 0.08534Como no hay interacción se sigue con la interpretacion de los dos factores. Al seguir con los factores los p.valores dieron mayor a 0.05, por tanto NO se rechaza la hipotesis nula y esto significa que fueron nulos los efectos de los dos factores (dosis del insecticida y numero de aplicaciones).
No es necesario hacer prueba de Tukey (prueba de comparacion de medias) ya que no hubo efecto ni de la dosis ni del numero de aplicaciones, lo que quiere decir que ni la dosis ni el numero de aplicaciones difieren. Esto estadisticamente hablando.
Revisión de supuestos
Prueba de normalidad de los residuales
\[Ho: \ Residuales \ normales \\ Ha: \ Residuales \ No \ normales\]
# Prueba de normalidad de los residuales
residual=mod9$residuals; residual## 2 10 16 4 13 6
## -0.48411735 -0.14761743 -0.45238058 -0.11607671 0.03207085 -0.25891959
## 17 8 14 5 9 1
## -0.09802343 -0.03590285 -0.05094965 0.41038930 0.10233013 0.23758720
## 12 11 15 7 3 18
## 0.21133006 0.25860398 -0.01651421 0.44641666 -0.07578647 0.03756007shapiro.test(residual)##
## Shapiro-Wilk normality test
##
## data: residual
## W = 0.96635, p-value = 0.727Los residuales son normales debido a que el p-valor dio mayor al 5%, por tanto se cumple el supuesto de normalidad. Esto tiene sentido porque los datos fueron generados con una distribucion normal.
Prueba de homocedasticidad de residuales
En cuanto a los diferentes tratamientos y sus varianzas \[Ho: Varianzas \ iguales \\ Ha: Al \ menos \ una \ diferente\]
# Prueba de homocedasticidad de residuales dosis
h=bartlett.test(residual~D$F1)
h##
## Bartlett test of homogeneity of variances
##
## data: residual by D$F1
## Bartlett's K-squared = 2.2328, df = 2, p-value = 0.3275El p-valor dio mayor al nivel de significancia del 5% por tanto, se cumple la igualdad de varianzas o supuesto de homocedasticidad.
# Prueba de homocedasticidad de residuales aplicaciones
bartlett.test(residual~D$F2)##
## Bartlett test of homogeneity of variances
##
## data: residual by D$F2
## Bartlett's K-squared = 0.39774, df = 2, p-value = 0.8197El p-valor dio mayor al nivel de significancia del 5% por tanto, se cumple la igualdad de varianzas o supuesto de homocedasticidad.
Supuesto de independencia de residuales
# Grafico de independencia de residuales
plot(residual, pch=16, col="purple", main="Gráfico de residuales")
Se observa que existe cierto patrón entre los residuales parecen estar relacionados, pero esto se puede deber debido a que la muestra es muy pequeña.
Gráfico de interacción
A continuación el gráfico de interacción de factores para comparar con los resultados del ANOVA:
library(ggplot2)
ggplot(data = D, aes(x = F1, y = biomasa, colour = F2,
group = F2)) +
stat_summary(fun = mean, geom = "point") +
stat_summary(fun = mean, geom = "line") +
labs(y = 'mean (D$biomasa)') +
theme_bw()
Como se puede observar en este gráfico presentado muestra que si existe interacción. Según el ANOVA dio que no había interacción, significa que estadisticamente no hay interaccion (guiandose por el p.valor). Para esta duda se va a seguir una propuesta que sugerió el profe en clase:
tabla=tapply(D$biomasa, list(D$F1, D$F2), mean); tabla## 4 5 6
## 3.25 2.965294 2.634592 2.592757
## 3.75 3.080186 3.128949 2.801410
## 4.25 3.344110 2.894908 3.021714X=addmargins(tabla, FUN=mean); X## Margins computed over dimensions
## in the following order:
## 1:
## 2:## 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.940436# con márgenes de la media# es lo mismo que arriba, acá solo están los márgenes
tabla1=tapply(D$biomasa, D$F1, mean); tabla1## 3.25 3.75 4.25
## 2.730881 3.003515 3.086911tabla2=tapply(D$biomasa, D$F2, mean); tabla2## 4 5 6
## 3.129864 2.886150 2.805294Acá se pueden observar las medias de las dosis y las del número de aplicaciones. Al mirar los márgenes se observa que la mejor dosis en cuanto al mayor rendimiento de biomasa sería una dosis de 4.25 (su promedio es el más alto 3.086911) y en cuanto a número de aplicaciones sería con 4 aplicaciones (su promedio sería 3.129864).De esto se concluiría que el mejor tratamiento sería una dosis de 4.25 con 4 aplicaciones. Ahora al momento de ver el cuerpo de la tabla, el mejor tratamiento sería una dosis de 4.25 con 4 aplicaciones (su promedio es el mayor con 3.344110). Al comparar los resultados se llega a lo mismo lo cual indica o es señal de que NO existe interacción entre los factores. Por tanto, es necesario hacer este paso cuando se tengan dudas o contradicciones con respecto al ANOVA y al GRÁFICO. En este caso el gráfico sirve para escoger cual fue el mejor tratamiento. Según este, el mejor tratamiento que se recomienda para obtener una mayor cantidad de biomasa es una dosis de 4.25 y aplicarla 4 veces durante el desarrollo del cultivo.
b. Presencia de covariable
set.seed(2020)
arcilla=sort(runif(18,0.20,0.40), decreasing= TRUE) # covariable
arcilla## [1] 0.3921445 0.3652331 0.3528828 0.3487672 0.3307115 0.3293806 0.3240412
## [8] 0.3237004 0.3079385 0.2953782 0.2845458 0.2818575 0.2788452 0.2786236
## [15] 0.2272194 0.2258305 0.2134769 0.2005165BIOMASA=D$biomasa; BIOMASA## [1] 2.708826 2.772692 2.143359 2.560519 2.708666 2.773705 2.770350 2.832470
## [9] 2.898280 3.359619 3.054099 3.157896 3.487669 3.451547 3.016111 3.042156
## [17] 3.200552 2.989329dfa=data.frame(D,BIOMASA,arcilla); dfa## F1 F2 biomasa BIOMASA arcilla
## 2 3.75 4 2.708826 2.708826 0.3921445
## 10 3.25 4 2.772692 2.772692 0.3652331
## 16 3.25 6 2.143359 2.143359 0.3528828
## 4 3.25 5 2.560519 2.560519 0.3487672
## 13 3.25 5 2.708666 2.708666 0.3307115
## 6 4.25 5 2.773705 2.773705 0.3293806
## 17 3.75 6 2.770350 2.770350 0.3240412
## 8 3.75 6 2.832470 2.832470 0.3237004
## 14 3.75 5 2.898280 2.898280 0.3079385
## 5 3.75 5 3.359619 3.359619 0.2953782
## 9 4.25 6 3.054099 3.054099 0.2845458
## 1 3.25 4 3.157896 3.157896 0.2818575
## 12 4.25 4 3.487669 3.487669 0.2788452
## 11 3.75 4 3.451547 3.451547 0.2786236
## 15 4.25 5 3.016111 3.016111 0.2272194
## 7 3.25 6 3.042156 3.042156 0.2258305
## 3 4.25 4 3.200552 3.200552 0.2134769
## 18 4.25 6 2.989329 2.989329 0.2005165#View(dfa)
plot(dfa$BIOMASA~dfa$arcilla,pch=16, main="Relación biomasa y arcilla", xlab="Arcilla", ylab="Biomasa (crecimiento)", col="red") # Grafico, se observa una relacion lineal entre covariable y variable respuesta (un supuesto)
reg_lin=lm(dfa$BIOMASA~dfa$arcilla)
abline(reg_lin)
Se observa que a medida que la cantidad de arcilla aumenta la cantidad de biomasa disminuye, parece que hay relación.
library(scatterplot3d)## Warning: package 'scatterplot3d' was built under R version 4.0.3scatterplot3d(dfa$arcilla,dfa$F1,BIOMASA, pch=16, color = "blue", type="h", main= "Relación biomasa, dosis y arcilla", zlab="Biomasa", xlab="Arcilla", ylab="Dosis")
Se observa que los puntos mas altos de biomasa se encuentran en las partes donde hay menor contenido de arcilla; en cuanto a la dosis parece ser que a medida que aumenta el crecimiento de biomasa es mayor. Lo que interesa en este ejercicio es saber si existe un efecto en la diminución de la biomasa a causa de la dosis de un insecticida. Se realizará el analisis respectivo.
Modelo del diseño factorial completo completamente al azar y en presencia de covariable
Meter el efecto de la covariable dentro del modelo Modelo de analisis de covarianza
\[y_{ijk} = \mu + \tau_i + \sigma_j + (\tau\sigma)_{ij} + \beta(x_{ijk}-\bar{x}) +\epsilon_{ijk}\\ i=1,2,3,\\ j=1,2,3;\\ k=1,2 \\ \ y \ sus \ condiciones \ laterales\]
Hipótesis para arcilla (covariable)
\[H_o:\beta\ = \ 0 \\ H_a:\beta\neq \ 0\]
modco=aov(BIOMASA~arcilla+F1*F2, data=dfa) #ANOVA
summary(modco)## 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 ' ' 1Acá tampoco hubo interacción entre los dos factores; en cuanto a F2 parece ser significativo al 6%. Lo que se puede observar es que se rechaza la hipotesis nula de la covariable arcilla debido que su p.valor dio menor a 0.05, de ahí que si hay efecto por parte de la arcilla.
modcos=aov(BIOMASA~arcilla+F1+F2, data=dfa) # Corriendo el modelo sin la interaccion
summary(modcos)## 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 ' ' 1Entonces si se justfica el uso de la covariable,lo que se puede observar es que el p.valor de la arcilla dio menor al nivel de significancia del 5%, de ahi que la cantidad de arcilla esta relacionada con la biomasa (se muestra un efecto significativo de la arcilla). En cuanto al F1 (dosis) al 5% No hay diferencias entre las dosis, pero si se encontró significancia en el F2 (número de aplicaciones). De ahí, que sí hay efecto por parte del número de aplicaciones.
El análisis de covarianza se entiende conociendo los fundamentos de la regresión simple y el análisis de varianza.
Fuente: https://www.datanovia.com/en/lessons/ancova-in-r/ Un análisis de covarianza permite comparar las medias ajustadas de dos o más grupos independientes.
SUPUESTOS DEL MODELO 1. Igualdad de pendientes: Debe existir relacion lineal entre la covariable x y la respuesta y, es decir \(\beta \neq 0\) (Linealidad entre la covariable y la variable de resultado)
summary(reg_lin) # para obtener el valor de la pendiente##
## Call:
## lm(formula = dfa$BIOMASA ~ dfa$arcilla)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.5982 -0.1471 -0.0296 0.1008 0.4786
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.0165 0.3699 10.858 8.64e-09 ***
## dfa$arcilla -3.6129 1.2228 -2.955 0.00932 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2748 on 16 degrees of freedom
## Multiple R-squared: 0.353, Adjusted R-squared: 0.3126
## F-statistic: 8.73 on 1 and 16 DF, p-value: 0.00932El valor de la pendiente o en el caso del modelo \(\beta\) es de -3.6129 , no es nula, lo que significa que hay una relacion lineal negativa entre la covariable y la variable respuesta biomasa. 2.Homogeneidad de las pendientes de regresión 3.Las covariables no son afectadas por los tratamientos 4.Los residuos provienen de una distribucion normal, o sea deben ser normales.(La variable de resultado debe tener una distribución aproximadamente normal)
qqnorm(resid(modcos), main="Normal Q-Q Plot")
qqline(resid(modcos), col="red")
Se observan los datos cerca a la línea, significa que los residuales son normales.
- Homoscedasticidad
- No hay valores atípicos significativos en los grupos.
Revisión de supuestos para el ANCOVA
\[Ho: \ Residuales \ normales \\ Ha: \ Residuales \ No \ normales\]
Prueba de normalidad de los residuales
# Prueba de normalidad de los residuales
resCo=modcos$residuals; resCo## 2 10 16 4 13 6
## -0.19373151 0.03959385 -0.24768346 -0.01999874 0.05102355 0.01868294
## 17 8 14 5 9 1
## -0.02830174 0.03236298 -0.14104993 0.26663743 0.27945907 0.06865770
## 12 11 15 7 3 18
## 0.29386790 0.06408276 -0.17529526 0.10840710 -0.27247069 -0.14424395shapiro.test(resCo)##
## Shapiro-Wilk normality test
##
## data: resCo
## W = 0.9418, p-value = 0.3109Los residuales son normales debido a que el p-valor dio mayor al 5%, por tanto se cumple el supuesto de normalidad.
Prueba de homocedasticidad de residuales
\[Ho: Varianzas \ iguales \\ Ha: Varianzas \ desiguales\]
# Prueba de homocedasticidad de residuales
bartlett.test(resCo, dfa$F1)##
## Bartlett test of homogeneity of variances
##
## data: resCo and dfa$F1
## Bartlett's K-squared = 1.8718, df = 2, p-value = 0.3922# Prueba de homocedasticidad de residuales
bartlett.test(resCo, dfa$F2)##
## Bartlett test of homogeneity of variances
##
## data: resCo and dfa$F2
## Bartlett's K-squared = 0.29936, df = 2, p-value = 0.861Los p-valor dieron mayor al nivel de significancia del 5% por tanto, se cumple la igualdad de varianzas o supuesto de homocedasticidad.
De las revisiones de supuestos presentadas anteriormente del ANOVA y del ANCOVA dieron que los residuales son normales, las varianzas son iguales y los residuales no son independientes. Lo que se percibe es que los valores de cada supuesto son parecidos.
No se recomienda el uso de arcilla debido a que cuando hay mayor cantidad es mas bajo el contenido o crecimiento de biomasa.
Gráfico de interacción en presencia de covariable
Este va a dar igual al primero porque para su representación se tienen en cuenta los factores y no la covariable
library(ggplot2)
ggplot(data = D, aes(x = F1, y = biomasa, colour = F2,
group = F2)) +
stat_summary(fun = mean, geom = "point") +
stat_summary(fun = mean, geom = "line") +
labs(y = 'mean (D$biomasa)') +
theme_bw()
EJERCICIO 5
- Existe un tipo de diseño anidado (factorial incompleta) 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. 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. Use los datos que se muestran para estimar las varianzas antes descritas. Una ayuda para la solución de este problema puede encontrarse en el libro: Design and Analysis of Experiments with R de John Lawson. ## a. Estimación de varianzas
library(readxl)
v <- read_excel("~/UNAL/R/PT.xlsx", sheet = "Hoja3")
View(v)
Data=data.frame(v); Data## respuesta ue finca variedad test
## 1 11.91 1 1 2 1
## 2 9.76 1 1 1 1
## 3 9.24 1 1 1 2
## 4 9.02 1 2 1 1
## 5 10.00 2 1 2 1
## 6 10.65 2 1 1 1
## 7 7.77 2 1 1 2
## 8 13.69 2 2 1 1
## 9 8.02 3 1 2 1
## 10 6.50 3 1 1 1
## 11 6.26 3 1 1 2
## 12 7.95 3 2 1 1
## 13 9.15 4 1 2 1
## 14 8.08 4 1 1 1
## 15 5.28 4 1 1 2
## 16 7.46 4 2 1 1
## 17 7.43 5 1 2 1
## 18 7.84 5 1 1 1
## 19 5.91 5 1 1 2
## 20 6.11 5 2 1 1
## 21 7.01 6 1 2 1
## 22 9.00 6 1 1 1
## 23 8.38 6 1 1 2
## 24 8.58 6 2 1 1
## 25 11.13 7 1 2 1
## 26 12.81 7 1 1 1
## 27 13.58 7 1 1 2
## 28 10.00 7 2 1 1
## 29 14.07 8 1 2 1
## 30 10.62 8 1 1 1
## 31 11.71 8 1 1 2
## 32 14.56 8 2 1 1
## 33 4.08 9 1 2 1
## 34 4.88 9 1 1 1
## 35 4.96 9 1 1 2
## 36 4.76 9 2 1 1
## 37 6.73 10 1 2 1
## 38 9.38 10 1 1 1
## 39 8.02 10 1 1 2
## 40 6.99 10 2 1 1
## 41 6.59 11 1 2 1
## 42 5.91 11 1 1 1
## 43 5.79 11 1 1 2
## 44 6.55 11 2 1 1
## 45 5.77 12 1 2 1
## 46 7.19 12 1 1 1
## 47 7.22 12 1 1 2
## 48 8.33 12 2 1 1
## 49 8.12 13 1 2 1
## 50 7.93 13 1 1 1
## 51 6.48 13 1 1 2
## 52 7.43 13 2 1 1
## 53 3.95 14 1 2 1
## 54 3.70 14 1 1 1
## 55 2.86 14 1 1 2
## 56 5.92 14 2 1 1
## 57 5.96 15 1 2 1
## 58 4.64 15 1 1 1
## 59 5.70 15 1 1 2
## 60 5.88 15 2 1 1
## 61 4.18 16 1 2 1
## 62 5.94 16 1 1 1
## 63 6.28 16 1 1 2
## 64 5.24 16 2 1 1
## 65 11.25 17 1 2 1
## 66 9.50 17 1 1 1
## 67 8.00 17 1 1 2
## 68 11.14 17 2 1 1
## 69 9.51 18 1 2 1
## 70 10.93 18 1 1 1
## 71 12.16 18 1 1 2
## 72 12.71 18 2 1 1
## 73 16.79 19 1 2 1
## 74 11.95 19 1 1 1
## 75 10.58 19 1 1 2
## 76 13.08 19 2 1 1
## 77 7.51 20 1 2 1
## 78 4.34 20 1 1 1
## 79 5.45 20 1 1 2
## 80 5.21 20 2 1 1library(lme4)## Loading required package: Matrix## Registered S3 methods overwritten by 'lme4':
## method from
## cooks.distance.influence.merMod car
## influence.merMod car
## dfbeta.influence.merMod car
## dfbetas.influence.merMod carData$respuesta<- as.numeric(as.character(Data$respuesta))
modelo <- lmer(respuesta ~ 1 +(1 | ue) + (1 | ue:finca) + +(1 | ue:finca:variedad), data = Data)## boundary (singular) fit: see ?isSingularsummary(modelo)## Linear mixed model fit by REML ['lmerMod']
## Formula: respuesta ~ 1 + (1 | ue) + (1 | ue:finca) + +(1 | ue:finca:variedad)
## Data: Data
##
## REML criterion at convergence: 326
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.92753 -0.39932 0.00922 0.43797 1.65397
##
## Random effects:
## Groups Name Variance Std.Dev.
## ue:finca:variedad (Intercept) 1.2305 1.1093
## ue:finca (Intercept) 0.0000 0.0000
## ue (Intercept) 7.0127 2.6481
## Residual 0.8795 0.9378
## Number of obs: 80, groups: ue:finca:variedad, 60; ue:finca, 40; ue, 20
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 8.2369 0.6188 13.31
## convergence code: 0
## boundary (singular) fit: see ?isSingularLa varianza atribuible a la parecela (ue) dio \(7.0127\) , la varianza atribuible a las parcelas anidadas en fincas dio \(0\) y la varianza de las variedades dentro de la finca dio \(1.2305\). De esto se puede concluir que existió mayor varianza en la parcela, por tanto esta es la que causa mayor variabilidad e imprecisión. Al aplicar la siguiente fórmula: \((100×(7.0127)/(7.0127+1.2305+0.8795))\) en donde \(0.8795\) se refiere a los residuales, se obtiene que el 76,8% de la variación total es debida a la variabilidad entre parcelas.
EJERCICIO 6.
- En el enlace https://cran.r-project.org/web/packages/asbio/asbio.pdf se tienen unos datos de potasio de muestras de suelos medidas en 8 diferentes laboratorios. Compare descriptivamente (medidas, tablas y gráficos) para representar los datos. ¿qué prueba me recomendaría para comparar la medida que usted seleccione? Proponga una solución. Sabiendo que son muestras mezcladas de una misma finca, ¿ Se perciben diferencias en las medidas como consecuencia probable de los laboratorios? Sugerencia: Use el enfoque no paramétrico considerado en clase y su respectiva prueba de comparación por pares (Nemenyi) ## Análisis descriptivo
library(asbio)## Warning: package 'asbio' was built under R version 4.0.3## Loading required package: tcltklibrary(lattice)
k=data(K)
df=data.frame(K)
#View(df)
potasio=df$K; potasio## [1] 296 260 341 359 323 321 287 413 335 315 330 326 354 266 348 343 284 324 351
## [20] 302 395 357 400 187 376 283 198 327 354 308 274 324 305 347 297 305 326 301
## [39] 316 312 297 280 300 319 286 218 280 241 226 243 199 205 225 227 338 303 341
## [58] 311 355 269 284 279 339 359 318 313 352 334 356 342 299 353laboratorios=df$lab; laboratorios## [1] B B B B B B B B B D D D D D D D D D E E E E E E E E E F F F F F F F F F G G
## [39] G G G G G G G H H H H H H H H H I I I I I I I I I J J J J J J J J J
## Levels: B D E F G H I Jplot(potasio, laboratorios, xlab="Contenido de K mg/kg de suelo", ylab="Laboratorios", pch=16, col = df$lab, main="Resultados muestras de K por laboratorio")
bwplot(potasio~laboratorios, ylab="Contenido de K mg/kg de suelo", xlab="Laboratorios", col=c(123,58,156,632,237,0,767,414), main="Resultados muestras de K por laboratorio")
# Ahora en valores el gráfico de caja
tapply(potasio, laboratorios, mean)## B D E F G H I J
## 326.1111 321.1111 316.5556 315.6667 304.1111 229.3333 313.2222 336.2222tapply(potasio, laboratorios, median)## B D E F G H I J
## 323 326 351 308 301 226 311 342tapply(potasio, laboratorios, var)## B D E F G H I J
## 1971.8611 855.8611 6456.2778 627.5000 236.3611 570.7500 985.6944 466.9444Los datos obtenidos permiten establecer que el laboratorio J fue el obtuvo mejores resultados en cuanto a la cantidad de K en el suelo, obtuvo una media de 336.2222 mg/kg suelo. En cuanto al laboratorio que obtuvo menores cantidades de K en el suelo fue el H con una media de 229.3333 mg/kg suelo. También se observa que los resultados obtenidos en el laboratorio E tuvieron una mayor variabilidad.
ANOVA
modk=aov(potasio~laboratorios); modk #busca la relacion entre las variables## Call:
## aov(formula = potasio ~ laboratorios)
##
## Terms:
## laboratorios Residuals
## Sum of Squares 68929.87 97370.00
## Deg. of Freedom 7 64
##
## Residual standard error: 39.00521
## Estimated effects may be unbalancedsmodk=summary(modk) #Resumen de la tabla del ANOVA
smodk## Df Sum Sq Mean Sq F value Pr(>F)
## laboratorios 7 68930 9847 6.472 8.92e-06 ***
## Residuals 64 97370 1521
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Al observar el p.valor de los laboratorios este dio menor al 5% de ahí que existen diferencias en los laboratorios.
Revisión de supuestos para el ANOVA
Prueba de normalidad de los residuales
\[Ho: \ Residuales \ normales \\ Ha: \ Residuales \ No \ normales\]
# Prueba de normalidad de los residuales
resk=modk$residuals; resk## 1 2 3 4 5 6
## -30.111111 -66.111111 14.888889 32.888889 -3.111111 -5.111111
## 7 8 9 10 11 12
## -39.111111 86.888889 8.888889 -6.111111 8.888889 4.888889
## 13 14 15 16 17 18
## 32.888889 -55.111111 26.888889 21.888889 -37.111111 2.888889
## 19 20 21 22 23 24
## 34.444444 -14.555556 78.444444 40.444444 83.444444 -129.555556
## 25 26 27 28 29 30
## 59.444444 -33.555556 -118.555556 11.333333 38.333333 -7.666667
## 31 32 33 34 35 36
## -41.666667 8.333333 -10.666667 31.333333 -18.666667 -10.666667
## 37 38 39 40 41 42
## 21.888889 -3.111111 11.888889 7.888889 -7.111111 -24.111111
## 43 44 45 46 47 48
## -4.111111 14.888889 -18.111111 -11.333333 50.666667 11.666667
## 49 50 51 52 53 54
## -3.333333 13.666667 -30.333333 -24.333333 -4.333333 -2.333333
## 55 56 57 58 59 60
## 24.777778 -10.222222 27.777778 -2.222222 41.777778 -44.222222
## 61 62 63 64 65 66
## -29.222222 -34.222222 25.777778 22.777778 -18.222222 -23.222222
## 67 68 69 70 71 72
## 15.777778 -2.222222 19.777778 5.777778 -37.222222 16.777778shapiro.test(resk)##
## Shapiro-Wilk normality test
##
## data: resk
## W = 0.94688, p-value = 0.004268Los residuales NO son normales debido a que el p-valor dio menor al 5%, por tanto NO se cumple el supuesto de normalidad.
Prueba de homocedasticidad de residuales
\[Ho: Varianzas \ iguales \\ Ha: Varianzas \ desiguales\]
# Prueba de homocedasticidad de residuales
bartlett.test(resk, laboratorios)##
## Bartlett test of homogeneity of variances
##
## data: resk and laboratorios
## Bartlett's K-squared = 32.201, df = 7, p-value = 3.727e-05E p-valor dio menor al nivel de significancia del 5% por tanto, NO se cumple la igualdad de varianzas o supuesto de homocedasticidad.
Análisis de varianza que permite varianzas desiguales Welch
library(car)
#Análisis de varianza que permite varianzas desiguales Welch
oneway.test(potasio~laboratorios, data=df, var.equal=FALSE)##
## One-way analysis of means (not assuming equal variances)
##
## data: potasio and laboratorios
## F = 14.464, num df = 7.000, denom df = 27.119, p-value = 1.114e-07Este resultado es similar al primer ANOVA en donde dio que el p.valor es 1.114e-07 menor al 0.05 de ahí que se rechaza que los tratamientos son iguales, por tanto son diferentes las medidas de cada tratamiento.
Prueba de independencia de residuales
# Grafico de independencia de residuales
plot(resk, pch=16, col= "pink", main="Gráfico de residuales")
Los residuales se comportan independiente.
PRUEBA DE NORMALIDAD Kruskall-Wallis
Como NO se obtuvo normalidad se va a aplicar la prueba de Kruskal-Wallis que es un método NO paramétrico.
\[H_o: \mu_{R_{var1}} = \mu_{R_{var2}} = \mu_{R_{var3}} = \mu_{R_{var4}} = \mu_{R_{var5}} = \mu_{R_{var6}} = \mu_{R_{var7}} = \mu_{R_{var8}} \\ H_a: H_o \ es\ falsa \ ,\ (al \ menos \ una \ es \ diferente)\]
kruskal.test(potasio~laboratorios)##
## Kruskal-Wallis rank sum test
##
## data: potasio by laboratorios
## Kruskal-Wallis chi-squared = 24.482, df = 7, p-value = 0.000937Como el p.valor dio pequeño eso significa que hay diferencias estadísticas en la variable respuesta para los distintos laboratorios; se rechaza la hipótesis nula de que los rangos promedios asignados a cada tratamiento son iguales.
Comparación de medias posterior a la prueba de Kruskal-Wallis. Se recomienda utilizar la prueba de pares de Nemenyi para ver si existen diferencias en los laboratorios.
PRUEBA DE NEMENYI
library(PMCMR)## Warning: package 'PMCMR' was built under R version 4.0.3## PMCMR is superseded by PMCMRplus and will be no longer maintained. You may wish to install PMCMRplus instead.posthoc.kruskal.nemenyi.test(potasio~laboratorios) # Prueba de Nemenyi## 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: potasio by laboratorios
##
## 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: noneEstadísticamente se reporta que en la mayoría no hubo diferencias entre los laboratorios al tener un p.valor mayor a 0.05. Hubo diferencias entre el laboratorio B y H, D y H, E y H, F y H y entre H y J. El laboratorio H es el que difiere con los demás. Como se pudo observar a pesar de que se trató de muestras mezcladas de una misma finca, se perciben diferencias estadísticas en las medidas como consecuencia de los laboratorios, las comparaciones con el laboratorio H. Entonces gracias a la prueba de Nemenyi se pudo reafirmar la solución a lo planteado.
Ejercicio 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.
Diseño de parcelas divididas en bloques completos
Para esto se usó un diseño (factorial completo) en parcelas divididas con bloques completos generalizados al azar y en presencia de covariable. Se evaluo el rendimiento del cultivo de plátano en un ensayo de campo en parcelas divididas con bloques completos y en presencia de una covariable, se utilizaron tres variedades de plátano (Colicero, Indio y Bocadillo), tres niveles de enmienda de Potasio y cuatro parcelas que actuaron como bloque. El diseño se puso en 4 bloques de tres parcelas principales, cada una de ellas dividida en 3 subparcelas. Las variedades fueron asignadas a las parcelas pricipales debido que para este estudio se hacía más fácil la aleatorización tomándolas como parcela principal, y por consiguiente los tratamientos de enmienda de Potasio a las subparcelas. Al colocar la enmienda como subparcela se tuvo como consecuencia mayores ventajas con respecto a las variedades al momento de correr el análisis de varianza. Adicionalmente se decidió colocar como covariable el pH en el suelo utilizado en cada unidad experimental. Y el objetivo que se tiene es estudiar los efectos de los factores, del bloqueo y de la covariable en cuestión.
Modelo estadístico con covariable y bloque
\[y_{i,j,k} = \mu + \alpha_{i} + \gamma_k + \tau_{k(i)} + \gamma_{j} + (\alpha\gamma)_{ij} + \beta(x_{ijk}-\bar{x})+\epsilon_{k(ij)} \\ i=1,2,3 \\ j=1,2,3 \\ k=1,2,3,4 \\ Con \ sus \ condiciones \ laterales\] ### Datos
library(readxl)
pd.p <- read_excel("~/UNAL/R/PT.xlsx", sheet = "Hoja6")
head(pd.p)## # A tibble: 6 x 5
## Observación Parcela Variedad Potasio Rendimiento
## <dbl> <chr> <chr> <dbl> <chr>
## 1 1 I Colicero 0 15.5
## 2 2 I Colicero 30 17.5
## 3 3 I Colicero 60 20.8
## 4 4 I Indio 0 20.5
## 5 5 I Indio 30 24.5
## 6 6 I Indio 60 30.2#View(pd.p)
set.seed(2021)
Parcela=as.factor(pd.p$Parcela); Parcela## [1] I I I I I I I I I II II II II II II II II II III
## [20] III III III III III III III III IV IV IV IV IV IV IV IV IV
## Levels: I II III IVVariedad=as.factor(pd.p$Variedad); Variedad## [1] Colicero Colicero Colicero Indio Indio Indio Bocadillo
## [8] Bocadillo Bocadillo Colicero Colicero Colicero Indio Indio
## [15] Indio Bocadillo Bocadillo Bocadillo Colicero Colicero Colicero
## [22] Indio Indio Indio Bocadillo Bocadillo Bocadillo Colicero
## [29] Colicero Colicero Indio Indio Indio Bocadillo Bocadillo
## [36] Bocadillo
## Levels: Bocadillo Colicero IndioPotasio=as.factor(pd.p$Potasio); Potasio## [1] 0 30 60 0 30 60 0 30 60 0 30 60 0 30 60 0 30 60 0 30 60 0 30 60 0
## [26] 30 60 0 30 60 0 30 60 0 30 60
## Levels: 0 30 60Rendimiento=as.numeric(pd.p$Rendimiento); Rendimiento## [1] 15.5 17.5 20.8 20.5 24.5 30.2 15.6 18.2 18.5 18.9 20.2 24.5 15.0 20.5 18.9
## [16] 16.0 15.8 18.3 12.9 14.5 13.5 20.2 18.5 25.4 15.9 20.5 22.5 12.9 13.5 18.5
## [31] 13.5 17.5 14.9 12.5 11.9 10.5pH=sort(runif(36,0.58,0.68), decreasing = TRUE) #covariable
df7=data.frame(pd.p, pH); df7## Observación Parcela Variedad Potasio Rendimiento pH
## 1 1 I Colicero 0 15.5 0.6782987
## 2 2 I Colicero 30 17.5 0.6758182
## 3 3 I Colicero 60 20.8 0.6755339
## 4 4 I Indio 0 20.5 0.6747738
## 5 5 I Indio 30 24.5 0.6745695
## 6 6 I Indio 60 30.2 0.6742709
## 7 7 I Bocadillo 0 15.6 0.6734292
## 8 8 I Bocadillo 30 18.2 0.6712144
## 9 9 I Bocadillo 60 18.5 0.6702091
## 10 10 II Colicero 0 18.9 0.6680486
## 11 11 II Colicero 30 20.2 0.6667538
## 12 12 II Colicero 60 24.5 0.6637490
## 13 13 II Indio 0 15.0 0.6620053
## 14 14 II Indio 30 20.5 0.6615422
## 15 15 II Indio 60 18.9 0.6614989
## 16 16 II Bocadillo 0 16.0 0.6583780
## 17 17 II Bocadillo 30 15.8 0.6509682
## 18 18 II Bocadillo 60 18.3 0.6501346
## 19 19 III Colicero 0 12.9 0.6440439
## 20 20 III Colicero 30 14.5 0.6436324
## 21 21 III Colicero 60 13.5 0.6403241
## 22 22 III Indio 0 20.2 0.6372207
## 23 23 III Indio 30 18.5 0.6367453
## 24 24 III Indio 60 25.4 0.6352266
## 25 25 III Bocadillo 0 15.9 0.6345698
## 26 26 III Bocadillo 30 20.5 0.6305494
## 27 27 III Bocadillo 60 22.5 0.6251267
## 28 28 IV Colicero 0 12.9 0.6192494
## 29 29 IV Colicero 30 13.5 0.6181744
## 30 30 IV Colicero 60 18.5 0.6068485
## 31 31 IV Indio 0 13.5 0.6066680
## 32 32 IV Indio 30 17.5 0.6051571
## 33 33 IV Indio 60 14.9 0.6024886
## 34 34 IV Bocadillo 0 12.5 0.5939581
## 35 35 IV Bocadillo 30 11.9 0.5832776
## 36 36 IV Bocadillo 60 10.5 0.5827267Esquema del diseño
library(collapsibleTree)
collapsibleTree(pd.p, hierarchy = c("Parcela","Variedad","Potasio"))Hipótesis
\[H_{01}:(\tau\beta)_{ij}=0;\forall_{ij}\]
\[H_{02}:\tau_i=0;\forall_{i}\] \[H_{03}:\beta_j=0;\forall_{j}\]
Análisis de varianza completo
library(lme4)
# Uso la función lmer
library(lmerTest)## Warning: package 'lmerTest' was built under R version 4.0.3##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## step#El rendimiento no es un factor
Potasios=as.factor(df7$Potasio); Potasios## [1] 0 30 60 0 30 60 0 30 60 0 30 60 0 30 60 0 30 60 0 30 60 0 30 60 0
## [26] 30 60 0 30 60 0 30 60 0 30 60
## Levels: 0 30 60Rendimientos=as.numeric(as.character(Rendimiento)); Rendimientos## [1] 15.5 17.5 20.8 20.5 24.5 30.2 15.6 18.2 18.5 18.9 20.2 24.5 15.0 20.5 18.9
## [16] 16.0 15.8 18.3 12.9 14.5 13.5 20.2 18.5 25.4 15.9 20.5 22.5 12.9 13.5 18.5
## [31] 13.5 17.5 14.9 12.5 11.9 10.5modelo.pdp7=lmer(Rendimientos~pH+Parcela+Variedad*Potasios+(1|Parcela:Variedad), data=df7)
anova(modelo.pdp7)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## pH 1.027 1.0274 1 18.4819 0.2216 0.643337
## Parcela 1.984 0.6614 3 7.1147 0.1427 0.931208
## Variedad 12.659 6.3296 2 7.4685 1.3652 0.312384
## Potasios 54.033 27.0167 2 20.3562 5.8270 0.009978 **
## Variedad:Potasios 8.005 2.0013 4 18.6550 0.4316 0.784009
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Como se puede observar, comenzando con la doble interacción entre variedad y potasio el p.valor (0.7676) dio mayor al nivel del 5% de significancia o 0.05 esto quiere decir que no hay interacción entre estas variables, de ahí que se puede seguir interpretando los p.valor. También se puede observar no tuvieron efecto los factores ni la covariable.
modelo.pdp7a=lmer(Rendimientos~pH+Parcela+Variedad+Potasios+(1|Parcela:Variedad), data=df7)
anova(modelo.pdp7a)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## pH 2.507 2.507 1 25.5235 0.6006 0.445478
## Parcela 2.017 0.672 3 7.4796 0.1610 0.919406
## Variedad 11.842 5.921 2 7.7768 1.4186 0.298385
## Potasios 67.157 33.578 2 23.8355 8.0445 0.002139 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Al correr el modelo sin interacción, se puede observar que solo la enmienda de potasio tuvo efecto, esto debido a que su p.valor dio menor a 0.06 (se esta trabajando con un nivel de significancia del 6%). Se observa que ni el pH, ni la parcela ni la variedad tuvieron efecto.
# sin la covariable
modelo.pdp7b=lmer(Rendimientos~Parcela+Variedad+Potasios+(1|Parcela:Variedad), data=df7)
anova(modelo.pdp7b)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Parcela 27.056 9.019 3 6 2.1836 0.1909357
## Variedad 12.880 6.440 2 6 1.5592 0.2849000
## Potasios 92.435 46.218 2 22 11.1901 0.0004442 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# sin la parcela
modelo.pdp7c=lmer(Rendimientos~Variedad+Potasios+(1|Parcela:Variedad), data=df7)
anova(modelo.pdp7c)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Variedad 9.236 4.618 2 9 1.1181 0.3683893
## Potasios 92.435 46.218 2 22 11.1901 0.0004442 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Como se observa el p.valor (0.0004442) del potasio dio significativo, menor a 0.05 de ahí que se rechaza la hipótesis nula de que su efecto es nulo como se comentó anteriormente, y por tanto se concluye que si hay efecto del potasio en cuanto al rendimiento, las enmiendas de potasio difieren.
Modelo estadístico sin covariable, sin bloque y sin interacción.
\[y_{i,j,k} = \mu + \alpha_{i} + \gamma_{j} +\epsilon_{ijk} \\ i=1,2,3 \\ j=1,2,3 \\ k=1,2,3,4 \\ Con \ sus \ condiciones \ laterales\] Factorial simple completamente al azar \[y_{i,j} = \mu + \alpha_{i} +\epsilon_{ij} \\ i=1,2,3 \\ j=1...12 \\ Con \ sus \ condiciones \ laterales\]
library(car)
# si se corre como factorial simple ya no hay efecto del potasio
MF1<- lm(Rendimientos ~ Potasios, data = df7); MF1##
## Call:
## lm(formula = Rendimientos ~ Potasios, data = df7)
##
## Coefficients:
## (Intercept) Potasios30 Potasios60
## 15.783 1.975 3.925anova(MF1)## Analysis of Variance Table
##
## Response: Rendimientos
## Df Sum Sq Mean Sq F value Pr(>F)
## Potasios 2 92.44 46.218 2.7957 0.07557 .
## Residuals 33 545.54 16.531
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1p.valor del potasio mayor a 0.06
Análisis descriptivo
library(lattice)
bwplot(Rendimientos~Potasios|Variedad, xlab = "Enmienda de Potasio kg/ha")
Como se puede observar en el gráfico al haber mayor cantidad de enmienda de potasio el rendimiento tiende a ser mayor en todas las variedades.El mejor tratamiento que se recomienda es el uso de la variedad Indio que como se observa es la que presenta más altos rendimientos y con una dosis o enmienda de potasio de 60 kg/ha, sin embargo este presenta mayor variabilidad. También se observó que el pH no tuvo efecto significativo, estos valores de pH se pueden considerar neutros, sin embargo se recomienda un pH neutro ya que estos funcionan para la mayoría de las variedades. Se sabe que el pH ótimo es 6.5. Por otro lado, es necesario ver los costos económicos que esto implica y la cantidad de dinero que tenga la persona al momento de cultivar.
EJERCICIO 8
Realice un resumen con la nota que aparece en las siguientes direcciones sobre: ## RESUMEN
El uso de los diseños en parcelas divididas
Un diseño de parcelas divididas se utiliza frecuentemente cuando se tiene un factor que es más difícil de aleatorizar que otros en campo. En ciencias relacionadas con la agronomía o el campo agrícola son de gran utilidad debido a que las parcelas divididas tuvieron origen en entornos agrícolas, su nombre viene de PARCELAS DE TIERRA. Sin embargo, este diseño también se ha visto utilizado en el entorno industrial. Una característica de nomenclatura de este diseño es que el factor de difícil aleatorización se denomina factor de parcela completa o parcela principal y los factores fáciles de aleatorizar se denominan factores de subparcela. El número total de tratamientos o experimentos corresponde al número de subparcelas, mientras que el número de parcelas completas corresponde al número de veces que se restablece los factores de la parcela principal. Ejemplos de factores difíciles de aleatorizar pueden ser: el cambio de temperatura de un horno, el tipo de riego en campo, la variedad de un cultivo, tipo de fertilizante entre otros. Esto también depende del investigador al momento de elegir como hacer el experimento y al momento de tener en cuenta los costos, ya que los presupuestos en algunos casos no son tan elevados entonces es necesario optimizar gastos. Ya existen paquetes analíticos en JMP y Design Expert que le ayudan al usuario al momento de elegir un buen experimento de parcela dividida y de hacer una interpretación adecuada de los datos. Es de gran importancia hacer un análisis cuidadoso debido a que se tienden a correr los datos (factores) como independientes cuando en realidad están correlacionados; esto haría que se concluya que el factor no es significativo cuando si lo era o al revés. Para finalizar, una ventaja de los experimentos de parcelas divididas es que tienden a producir estimaciones de parámetros del modelo precisas que le pueden ser de gran utilidad al investigador y que también le brindan mayor cantidad de información sobre su experimento.
Sobre lo que significa unidad experimental y unidad de observación
La unidad experimental también se conoce como unidad de replicación y es el ente más pequeño que se asigna independientemente de las demás unidades a un tratamiento particular. Se puede resumir como el ente que aporta el dato y que se le aplica un tratamiento de forma independiente, de ahí que es aquella que recibe el tratamiento algunos ejemplos de unidades experimentales pueden ser las vacas, plantas, parcelas, hojas, frutos, entre otros. Por otro lado, existe el término unidad observatorio o unidad de muestreo y este es un ente físico sobre el cual se mide un resultado de interés en un experimento; cuando se trata de un diseño simple estos dos términos mencionados anteriormente suelen ser iguales (mismo ente). El ente puede desempeñar el papel de los dos dentro de un mismo experimento (en parcelas divididas). Una estructura de datos jerárquica se refiere al arreglo de datos donde las observaciones no son independientes entre sí, sino que existe una correlación entre ellas. Por otro lado, cuando las unidades de observación se aniden dentro de una unidad experimental, estas se denominarán submuestras lo que indica que las observaciones están correlacionadas. La replicación involucra una repetición independiente de un tratamiento. Importancia de identificar las variables independientes que afectan a la variable respuesta, especificar la restricción en la aleatorización. Uso de modelos mixtos que reconocen a la vez fuentes de variabilidad aleatoria en un conjunto de datos. Para obtener un ANOVA correcto los valores promedio de las unidades de muestreo deben calcularse para cada unidad experimental. Si no se reconoce bien cuál es la unidad de muestreo se genera la pseudo replicación en donde se trata cada unidad de muestreo como si fuera una unidad experimental (se llega a una mala interpretación de los resultados).
Guía para diseñar experimentos exitosos
Es necesario al momento de diseñar un experimento tener varios pasos presentes, entre ellos está el ciclo de retroalimentación en donde le permite al investigador crear nuevas hipótesis científicas y modificaciones del diseño para eventos futuros. El diseño experimental hace referencia al conjunto de procedimientos y reglas por los cuales los tratamientos son asignados a las unidades experimentales. Los conceptos claves del diseño experimental son replicación (aplicar cada tratamiento a unidades experimentales múltiples y mutuamente independientes, obtención de variabilidad), aleatorización (asignar tratamientos a u.e de manera que cada una tenga la misma probabilidad de recibir cada tratamiento), bloqueo (grupo de u.e homogéneos), unidades experimentales (ente más pequeño al que se le aplica un tratamiento), unidad de observación (ente sobre el que se realizan las mediciones u observaciones), factor(tipo de tratamiento), nivel (estados de un factor).
REPLICACIÓN: permite estimar el error experimental, sin replicación es imposible estimar variación en las estimaciones de los efectos del tratamiento; con este se aumenta la precisión de un experimento; se le aplica al nivel de la unidad experimental. Existe la pseudorreplicación, en esta los tratamientos experimentales se replican, pero a una escala insuficiente para una interpretación válida de los resultados.
Primero definir la unidad experimental, hacer repeticiones su número tiene un efecto directo en el experimento es el tema más importante y estas pueden darse en cuatro niveles dentro del experimento: en la unidad experimental, en el experimento completo, en muestreo dentro de las unidades experimentales o en medidas repetidas. Es necesario identificar el tipo de distribución de la variable de interés; elegir una correcta estructura del diseño experimental que contenga bloques y varios niveles de replicación.
Diseños aumentados implican mayor cantidad de tratamientos no repetidos organizados en bloques pequeños. Se dice que los diseños más simples son los menos robustos, cada vez que se aumenta su grado de complejidad el diseño se vuelve más eficiente.
Los tratamientos que tienen distancias medias altas entre réplicas tenderán a tener varianzas más altas que con los que tienen distancias medias bajas entre réplicas (van Es y van Es, 1993).
Aleatorización en un experimento significa que es igualmente probable que cada tratamiento se aplique a cada unidad experimental y su uso principal es para minimizar el efecto de otras variables que no son las de interés., sirve como un seguro contra perturbaciones que se puedan generar dentro del experimento.
Bloqueo se utiliza por precisión o por conveniencia, un ejemplo que cabe resaltar es el de parcelas divididas que se utiliza cuando existe un factor que es difícil de aleatorizar ya sea porque requiere de unidades experimentales más grandes (ejemplos: tratamientos de riego, de labranza y fechas de siembra) por consiguiente se coloca en la parcela principal. Existen diseños en bloque incompletos como: cuadrado de celosia, latino; los diseños alfa, en columnas y en hileras ofrecen flexibilidad para variar el número de repeticiones, de tratamientos y el tamaño del bloque. Existen diseños en bloques incompletos y parcialmente balanceados. Tamaño de las unidades experimentales: entre más grande se logra una disminución en la varianza.
Es necesario poco a poco ir practicando y diseñando poco a poco experimentos y ahondar más sobre el tema si uno quiere lograr excelentes experimentos. El fracaso no es una opción en diseño de experimentos debido a que se requiere invertir un gran presupuesto en estos de ahí que hay que hacer lo mejor posible para pasar la prueba y no sentir que la plata se perdió.
REFERENCIAS http://www.ru.ac.bd/stat/wp-content/uploads/sites/25/2019/03/502_07_00_Lawson_Design-and-Analysis-of-Experiments-with-R-2017.pdf FloSchomo.2018.Error in REML code of running linear regression model equation. URL: https://stackoverflow.com/questions/53236333/error-in-reml-code-of-running-linear-regression-model-equation http://www.ru.ac.bd/stat/wp-content/uploads/sites/25/2019/03/502_07_00_Lawson_Design-and-Analysis-of-Experiments-with-R-2017.pdf FloSchomo.2018.Error in REML code of running linear regression model equation. URL: https://stackoverflow.com/questions/53236333/error-in-reml-code-of-running-linear-regression-model-equation