Parcial_DE_2020_2
Marly Daniela Quiroga Buitrago- 1056505936
8/11/2020
PUNTO 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 863. Con esta información se completó la tabla del ANOVA que se muestra a continuación:
Para completar la tabla anterior se tuvo en cuenta lo siguiente:
- En primer lugar se establecieron los grados de libertad \(Df\), para los tratamientos están dados por la ecuación \((a-1)\), donde \(a\) corresponde al número de tratamientos, como son 9 entonces los grados de libertad entre tratamientos son 8. Por otra parte, los grados de libertad \(Df\) del error están dados por \((N-a)\), donde \(N\) corresponde al número de observaciones totales (72) y \(a=9\), entonces el error presenta 63 grados de libertad.
a=9
N=72- Después se calculó la suma de cuadrados totales \(SC_{Totales}\) partiendo de que \(S^2=\frac{SC_{Totales}}{Df_{Totales}}\), como la varianza es \(S^2=863\) y los \(Df_{Totales}=71\), se despeja \(SC_{Totales}=S^2.{Df_{Totales}}\), entonces \(SC_{Totales}=61273\).
Df_Totales=71
Varianza=863
SC_T=(Varianza*Df_Totales)
SC_T## [1] 61273- A partir de \(SC_{Totales}=61273\) restándole la suma de cuadrados entre tratamientos \(SC_{Trat}=6000\) se calcula la suma de cuadrados del error \(SC_{Error}=55273\).
SC_Trat=6000
SC_Error=SC_T-SC_Trat
SC_Error## [1] 55273- Los cuadrados medios entre tratamientos se da por \(CM_{trat}=\frac{SC_{Trat}}{(a-1)}=750\) y los cuadrados medios del error es \(CM_{Error}=\frac{SC_{Error}}{(N-a)}=877.3492\).
CM_trat=SC_Trat/(a-1)
CM_trat## [1] 750CM_Error=SC_Error/(N-a)
CM_Error## [1] 877.3492- Finalmente, el F calculado está dado por \(\frac{CM_{trat}}{CM_{Error}}=0.855\)
F_calculado=CM_trat/CM_Error
F_calculado## [1] 0.8548478Se plantea la hipótesis nula: Igualdad de promedios de los indices de clorofila en todos los tratamientos:
\[H_0=\mu_{Trt1}=\mu_{Trt2}=\cdots=\mu_{Trt8}=\mu_{Trt9}\]
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 así como el cociente F calculado de la tabla para concluir)?
- Como no tenemos p-valor, una forma de calcularlo a partir de nuestro F calculado es con la función pf, como se muestra:
p_valor= pf(F_calculado,8,63, lower.tail = F)
p_valor## [1] 0.5588979- Teniendo en cuenta que el \(F_{Calculado}<F_{Tabulado}=0.855<2.8\), se acepta la hipotesis nula \(H_0\), esto quiere decir que el promedio del indice de clorofila en todos los tratamientos es igual. Lo anteior se confirma con el p-valor (0.559>0.05).
Curva F
F_tabulado<- 2.8
x <- seq( -4, 4, by = 0.1)
y <- dnorm( x )
plot( function(x) df( x, df1 = 8, df2= 63), 0, 4, ylim = c( 0, 1 ),
col = "black", type = "l", lwd = 2,
main = "Curva F" )
abline(v=F_tabulado,col="blue")
text(3.1,0.1,"f_tab=2.8")
abline(v=F_calculado, col = "red")
text(1.2,0.1,"f_cal=0.855")
text(3.5, 0.5, "Rechazo Ho")
text(2, 0.5, "No rechazo Ho")
¿vale la pena comparar las medias de tratamientos a posteriori del ANOVA (prueba de Tukey)?
En este caso NO es necesario aplicar esta prueba, dado que no se rechazó la hipótesis nula, indicando que las medias de los tratamientos son iguales. Lo anterior teniendo en cuenta que ya que la prueba de tukey se aplica cuando se rechaza la hipótesis de igualdad de medias en la prueba Anova.
PUNTO 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. Los datos se muestran en la siguiente tabla:
library(readxl)
Algodon <- read_excel("C:/Users/CK0001/Desktop/ING.AGRONOMICA.UNAL/semestre 6/D.E/Parcial/Algodon.xlsx") # Se organizó de otra manera la tabla en un archivo de excel
Algodon## # A tibble: 30 x 3
## Metodo Granjero perdidas
## <chr> <dbl> <dbl>
## 1 M1 1 NA
## 2 M1 2 6.75
## 3 M1 3 13.0
## 4 M1 4 10.3
## 5 M1 5 8.01
## 6 M1 6 8.42
## 7 M2 1 5.54
## 8 M2 2 3.53
## 9 M2 3 11.2
## 10 M2 4 7.21
## # ... with 20 more rowsAhora se convierte el archivo de excel en un data-frame y tanto los metodos como los granjeros se pasan a factores y las pérdidas a un vector, como se muestra a continuación:
Algodon_DQ = as.data.frame(Algodon)
Metodo_DQ = factor(Algodon$Metodo)
Granjero_DQ = factor(Algodon$Granjero)
Perdidas_DQ = as.vector(Algodon$perdidas)
str(Algodon_DQ)## 'data.frame': 30 obs. of 3 variables:
## $ Metodo : chr "M1" "M1" "M1" "M1" ...
## $ Granjero: num 1 2 3 4 5 6 1 2 3 4 ...
## $ perdidas: num NA 6.75 13.05 10.26 8.01 ...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)
Como es un modelo desbalanceado, para hacer el analisis de varianza ANOVA, se usa la función lm ya que esta corre este tipo de datos, el desbalance se da porque las mediciones de la granja 1 procesada por el limpiador M1 se perdieron. En este caso se usa en ANOVA tipo III, ya que uno de sus pros es que este no depende del tamaño de la muestra, es decir, es apto para los datos no balanceados donde tenemos un número desigual de observaciones en cada grupo y este es nuestro caso ya que uno de los datos se perdió. Entonces el ANOVA es:
library(rstatix)## Warning: package 'rstatix' was built under R version 4.0.3##
## Attaching package: 'rstatix'## The following object is masked from 'package:stats':
##
## filterModelo_DQ1 = lm(Perdidas_DQ ~ Granjero_DQ + Metodo_DQ) # Modelo lineal para modelos desbalanceados
ANOVA_DQ1 = Anova(Modelo_DQ1, type = "3") # Desarrollo del ANOVA desbalanceado mediante el uso de anova type = 3.
ANOVA_DQ1## Anova Table (Type III tests)
##
## Response: Perdidas_DQ
## Sum Sq Df F value Pr(>F)
## (Intercept) 198.62 1 143.8751 2.610e-10 ***
## Granjero_DQ 139.66 5 20.2331 5.163e-07 ***
## Metodo_DQ 49.12 4 8.8951 0.0003186 ***
## Residuals 26.23 19
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Modelo_DQ2 = lm(Perdidas_DQ ~ Metodo_DQ+Granjero_DQ) # Modelo lineal para modelos desbalanceados
ANOVA_DQ2 = Anova(Modelo_DQ2, type = "3") # Desarrollo del ANOVA desbalanceado mediante el uso de anova type = 3.
ANOVA_DQ2## Anova Table (Type III tests)
##
## Response: Perdidas_DQ
## Sum Sq Df F value Pr(>F)
## (Intercept) 198.62 1 143.8751 2.610e-10 ***
## Metodo_DQ 49.12 4 8.8951 0.0003186 ***
## Granjero_DQ 139.66 5 20.2331 5.163e-07 ***
## Residuals 26.23 19
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Al modificar el orden de colocación de los efectos del modelo en el ANOVA no afecta el resultado final, solo se observa que en la tabla resultante hay cambio de posición, es decir, en el ANOVA_DQ1 aparce primero Granjero y luego el Método, mientras en el ANOVA_DQ2 ocurre lo contrario. Sin embargo, lo que realmente importa que son los resultados del análisis no cambian.
Por otra parte, de acuerdo a los resultados de estos ANOVAS, se observa que existen diferencias significativas entre los metodos empleados para la limpieza del algodón y a su vez hay diferencias significativas entre las 6 granjas que lo producen.
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 limpiador de las fibras de algodón. Compare este resultado con el caso desbalanceado (de ser posible).
El dato que se perdió correspondiente al Granjero 1 con el método 1, se reemplazará con el promedio obtenido de los datos de los 5 granjeros con el mismo metodo (M1), para realizar el nuevo análisis de varianza para este modelo que ya no es desbalanceado, se crea nuevamente el archivo de excel con los datos completos, los datos son:
library(readxl)
Algodon2 <- read_excel("C:/Users/CK0001/Desktop/ING.AGRONOMICA.UNAL/semestre 6/D.E/Parcial/Algodon2.xlsx") # Se sube la tabla con todos los datos
Algodon2## # A tibble: 30 x 3
## Metodo Granjero perdidas
## <chr> <dbl> <dbl>
## 1 M1 1 9.30
## 2 M1 2 6.75
## 3 M1 3 13.0
## 4 M1 4 10.3
## 5 M1 5 8.01
## 6 M1 6 8.42
## 7 M2 1 5.54
## 8 M2 2 3.53
## 9 M2 3 11.2
## 10 M2 4 7.21
## # ... with 20 more rowsAhora se convierte el archivo de excel en un data-frame y tanto los metodos como los granjeros se pasan a factores y las pérdidas a un vector, como se muestra a continuación:
Algodon_DQ2 = as.data.frame(Algodon2)
Metodo_DQ2 = factor(Algodon2$Metodo)
Granjero_DQ2 = factor(Algodon2$Granjero)
Perdidas_DQ2 = as.vector(Algodon2$perdidas)
str(Algodon_DQ2)## 'data.frame': 30 obs. of 3 variables:
## $ Metodo : chr "M1" "M1" "M1" "M1" ...
## $ Granjero: num 1 2 3 4 5 6 1 2 3 4 ...
## $ perdidas: num 9.3 6.75 13.05 10.26 8.01 ...Se plantea la hipótesis: Las pérdidas medias de peso para los cinco métodos de limpiador de las fibras de algodón son iguales.
\[H_0=\mu_{Pérdidas M1}=\mu_{Pérdidas M2}=\mu_{Pérdidas M3}=\mu_{Pérdidas M4}=\mu_{Pérdidas M5}\]
Nuevo análisis de varianza para el modelo balanceado.
Modelo_DQ3 = lm(Perdidas_DQ2 ~ Granjero_DQ2+Metodo_DQ2) # Modelo lineal.
ANOVA_DQ3 = aov(Modelo_DQ3) # Desarrollo del ANOVA desbalanceado mediante el uso de anova type = 3.
summary(ANOVA_DQ3)## Df Sum Sq Mean Sq F value Pr(>F)
## Granjero_DQ2 5 139.36 27.873 21.015 2.45e-07 ***
## Metodo_DQ2 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 ' ' 1Sabiendo que el orden de los factores no afecta el resultado final, no es necesario hacer el ANOVA con el otro cambiando de lugar los factores. En esta tabla de ANOVA, se observa que también hay diferencias significativas entre los granjeros y los métodos usados en la limpieza del algodón. Comparando este análisis con el del modelo desbalanceado donde también hay diferencias significativas, entonces se concluye que con los dos análisis de ANOVA se rechaza la hipóteis nula que indica que las pérdidas medias de peso para cada granja y los cinco métodos de limpiador de las fibras de algodón son las mismas.
Para intentar de identificar cual método es mejor, es decir, el que ocasiona menor pérdida de peso en el algodón, se realizará un análisis descriptivo.
-Gráfico de cajas
library(ggplot2)
ggplot(Algodon_DQ2, aes(Metodo_DQ2, y = Perdidas_DQ2))+
geom_boxplot()
De acuerdo al gráfico uno de los mejores métodos es el M2, ya que este es el que ocasiona menos perdidas en el algodón al momento de limpiarlo.
Ahora visualizaremos la tabla de medias, para comparar el efecto de los métodos de limpieza sobre el peso del algodón respecto a las fincas:
M_perdidas <- tapply(Perdidas_DQ2, list(Metodo_DQ2, Granjero_DQ2), mean)
addmargins(M_perdidas, FUN= mean)## Margins computed over dimensions
## in the following order:
## 1:
## 2:## 1 2 3 4 5 6 mean
## M1 9.2980 6.750 13.05 10.260 8.010 8.420 9.298000
## M2 5.5400 3.530 11.20 7.210 3.240 6.450 6.195000
## M3 7.6700 4.150 9.79 8.270 6.750 5.500 7.021667
## M4 7.8900 1.970 8.97 6.120 4.220 7.840 6.168333
## M5 9.2700 4.390 13.44 9.130 9.200 7.130 8.760000
## mean 7.9336 4.158 11.29 8.198 6.284 7.068 7.488600De acuerdo a la tabla de medias, el mejor tratamiento (el que presenta menos perdidas) es el metodo M4 con el granjero 2, sin embargo, observando los margenes el mejor es el M2 con el granjero 2. Para concluir con esta situcación, haré un analisis de medias y respectivo coeficiente de variación:
medias_perdidas <- tapply(Perdidas_DQ2, list(Metodo_DQ2), FUN = mean)
medias_perdidas## M1 M2 M3 M4 M5
## 9.298000 6.195000 7.021667 6.168333 8.760000desviacion_perdidas=tapply(Perdidas_DQ2, list(Metodo_DQ2), sd)
coef_variacion = 100 * desviacion_perdidas/medias_perdidas
coef_variacion## M1 M2 M3 M4 M5
## 23.52622 47.01110 28.68060 42.91226 33.95876Hay que tener cuidado con los coeficientes de variación ya que son mayores al 20%, y esto puede ser peligroso para tomar decisiones, por lo que se recomendaría reevaluar el diseño. Pero haciendo análisis con estos CV, se observa que el M2 y M4 son los que presentan mayor variación, entonces por esta razón se descartan inmediatamente, entonces, el M1 presenta el menor CV pero también es el que ocasiona mayores perdidas de peso en el algodon (descartado). Finalmente, el M3 es el que sigue en menores perdidas de acuerdo a la tabla de medias y a su vez presenta un menor CV, por lo que este podría ser el mejor método.
PUNTO 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.6+0.1 y de la capa inferior osciló entre 2 y 2.5+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. Use la función sort.int de R para ordenar los datos de cada capa con la opción partial=25+6 dentro de la propia función sort.int.
set.seed(2021)
CO_5 = runif(50,min=3,max=3.7) # Generación de datos uniformes para la capa superior
CO_10 = runif(50,min=2,max=2.7) # Para la capa inferior
CO5=sort.int(CO_5, partial = 31)
CO10=sort.int(CO_10, partial = 31)- se crea la variable ventana_DQ y se definen los puntos para generar los 50 datos
Ventana_DQ=expand.grid(longitud= seq(0,100,25), latitud = seq(0,200,length.out = 10))
plot(Ventana_DQ)
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.
- Ahora se crea el data.frame para el gráfico 3D
DQ <- data.frame(longitud = rep(Ventana_DQ$longitud,2),
latitud = rep(Ventana_DQ$latitud,2),
profundidad = rep(c(-5,-10),each = 50),
col= c(CO5,CO10))- Se hace el gráfico 3D con la función plot_ly:
library(plotly)## Warning: package 'plotly' was built under R version 4.0.3##
## 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=DQ$longitud, y=DQ$latitud, z=DQ$profundidad, type="scatter3d", mode="markers", color=DQ$co)## 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.Compare si se encuentran diferencias en la media de carbono entre capas utilizando un nivel de confianza del 95%.
Como los niveles de carbono de la capa superior influye sobre los niveles de la capa inferior (están relacionados) se realiza una prueba t-pareada:
prueba_co= t.test(DQ$co~DQ$profundidad, alternative = "g", paired = T); prueba_co##
## Paired t-test
##
## data: DQ$co by DQ$profundidad
## t = -42.878, df = 49, p-value = 1
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
## -1.121879 Inf
## sample estimates:
## mean of the differences
## -1.079664ifelse(prueba_co$p.value<0.05,"Rechazo Ho","No rechazo Ho")## [1] "No rechazo Ho"- Como el p-valor 2.2e-16<0.05, si se encuentran diferencias en la media de carbono entre las capas utilizando un nivel de confianza del 95%. De acuerdo al gráfico, a mayor profundidad, los valores de CO disminuyen tal cual ocurre en la realidad.
PUNTO 4. Diseño factorial completo \((3^2)\) en arreglo completamente al azar
Este diseño está conformado por 2 factores (Dosis del insecticida y Número de aplicaciones), cada uno con 3 niveles, es decir, presenta 9 combinaciones de tratamientos cada uno se repite dos veces para un total de 18 observaciones. Los factores y la respuesta fueron creados con el código:
D <- expand.grid(Dosis_insecticida=c(3.25, 3.75, 4.25), Numero_aplicaciones=c(4,5,6), Repeticiones=paste0("r",1:6)) #Crea el diseño 3^2
D <- rbind(D, D) # Crea la estructura para dos repeticiones por tratamiento
set.seed(2020)
D <- D[order(sample(1:18)),] # Aleatoriza la estructura
class(D)## [1] "data.frame"D$biomasa=sort.int(rnorm(18,3,0.3),partial = 9) #Crea la respuesta
D## Dosis_insecticida Numero_aplicaciones Repeticiones biomasa
## 2 3.75 4 r1 2.708826
## 10 3.25 4 r2 2.772692
## 16 3.25 6 r2 2.143359
## 4 3.25 5 r1 2.560519
## 13 3.25 5 r2 2.708666
## 6 4.25 5 r1 2.773705
## 17 3.75 6 r2 2.770350
## 8 3.75 6 r1 2.832470
## 14 3.75 5 r2 2.898280
## 5 3.75 5 r1 3.359619
## 9 4.25 6 r1 3.054099
## 1 3.25 4 r1 3.157896
## 12 4.25 4 r2 3.487669
## 11 3.75 4 r2 3.451547
## 15 4.25 5 r2 3.016111
## 7 3.25 6 r1 3.042156
## 3 4.25 4 r1 3.200552
## 18 4.25 6 r2 2.989329Ahora se convierten las variables a factores y la respuesta a vector, seguido de la creación de un data.frame con estos datos como se muestra a continuación:
Dosis_insecticida1=as.factor(D$Dosis_insecticida)
Numero_aplicaciones1 = as.factor(D$Numero_aplicaciones)
Biomasa1 = as.vector(D$biomasa)
D1<-data.frame(Dosis_insecticida1,Numero_aplicaciones1,Biomasa1)
str(D1)## 'data.frame': 18 obs. of 3 variables:
## $ Dosis_insecticida1 : Factor w/ 3 levels "3.25","3.75",..: 2 1 1 1 1 3 2 2 2 2 ...
## $ Numero_aplicaciones1: Factor w/ 3 levels "4","5","6": 1 1 3 2 2 2 3 3 2 2 ...
## $ Biomasa1 : num 2.71 2.77 2.14 2.56 2.71 ...Por otra parte, se realiza un pequeño análisis descriptivo para tener una idea de como es el comportamiento de los datos e inferir maso menos con cual tratamiento se afecta menos el crecimiento del cultivo, es decir, con cual tratamiento el cultivo presenta una mayor biomasa, entonces:
lattice::bwplot(D1$Biomasa1~D1$Numero_aplicaciones1|D1$Dosis_insecticida1)
De acuerdo al gráfico se puede interpretar que el tratamiento que tiene menor efecto en el crecimiento del cultivo, es decir, donde las plantas presentan mayor biomasa, se alcanza en la dosis de 4.25 con 4 aplicaciones, además, se observa que existe poca variabilidad entre estos datos respecto a otros tratamientos, por ahora se podría considerar el tratamiento 4.25 con 4 aplicaciones como una buena opción para el manejo de malezas sin que se afecte en gran medida el cultivo.
Una forma para visualizar el diseño es:
library(collapsibleTree) ## Warning: package 'collapsibleTree' was built under R version 4.0.3collapsibleTree(D,hierarchy=c("Dosis_insecticida", "Numero_aplicaciones","Repeticiones"))a. Escriba (completamente especificado) el Modelo del diseño
\[y_{ijk}=\mu+\alpha_i+\beta_j+(\alpha\beta)_{ij}+\epsilon_{ijk}\\i:1\cdots3\\j:1,3\\k:1\cdots6\]
Donde:
\(y_{ijk}\) Es la \(ikj-esima\) biomasa para el \(i-esimo\) nivel del factor 1 (Dosis_insecticida) y el \(j-esimo\) nivel del factor 2 (Numero_aplicaciones).
\(\mu\) Es la media general de Biomasa
\(\alpha_i\) Es el efecto del \(i-esimo\) nivel del factor 1 (Dosis_insecticida).
\(\beta_j\) Es el efecto del \(j-esimo\) nivel del factor 2 (Numero_aplicaciones).
\((\alpha\beta)_{ij}\) Es la interacción del \(i-esimo\) nivel del factor 1 (Dosis_insecticida) con el \(j-esimo\) nivel del factor 2 (Numero_aplicaciones).
\(\epsilon_{ijk}\) Es el error residual para el \(i-esimo\) nivel del factor 1 (Dosis_insecticida), \(j-esimo\) nivel del factor 2 (Numero_aplicaciones) y la \(k-esima\) repetición.
Plantear Hipótesis
El efecto de la interacción entre los dos factores es nulo para todo i,j. \[H_{01}:(\alpha\beta)_{ij}=0~\\ \forall i.j\]
El efecto del F1, en este caso el efcto de la la Dosis del insecticida es nulo. \[H_{02}:\alpha_i=0~\\ i=1\cdots3\]
El efecto del F2, en este caso el efecto del número de aplicaciones es nulo. \[H_{03}: \beta_j=0~\\ j=1,3\]
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.
ANOVA
# Calculo de la tabla Anova
modelo1 <- aov(D1$Biomasa1~(D1$Dosis_insecticida1*D1$Numero_aplicaciones1))
summary(modelo1)## Df Sum Sq Mean Sq F value Pr(>F)
## D1$Dosis_insecticida1 2 0.4161 0.20804 1.979 0.194
## D1$Numero_aplicaciones1 2 0.3426 0.17128 1.630 0.249
## D1$Dosis_insecticida1:D1$Numero_aplicaciones1 4 0.1635 0.04087 0.389 0.812
## Residuals 9 0.9459 0.10510De acuerdo alanálisis de varianza, no existe interacción entre las dosis del insecticida y el numero de aplicaciones, dado que se observa que segun el pvalor (0.812>0,05), entoceS, segun lo anterior se acepta la primer hipotesis nula \(H_{01}\), la cual dice que la interacción entre los dos factores es nula. Por otra parte, observando el F-value, este presenta valores mayores que 1 en los factores, lo que indica que la principal fuente de variación son los tratamientos y no las replicaciones, lo cual se supone es lo deseable e ideal, porque si se presentaran valores menores que 1 es preocupante para el investigador, ya que indicaria que la principal fuentes de variación son las replicasciones y esto no debería ser así, entonces en ese caso debería reevaluarse el diseño.
En cuanto a la segunda y tercera hipótesis, que indican que el efecto de los dos factores (Dosis y número deaplicaciones) es nulo sobre la biomasa, en este caso se observan p-valores para el factor Dosis de 0.194 y para el factor numero de palicaciones de 0.249, los cuales son mayores de 0.05, quiere decir que los resultados no son significativos, entonces tampoco se rechazan las hipóteis por lo que el efecto del Factor 1 y del Factor 2 son nulos.
Prueba de comparación de medias: Se puede realizar debido a que según el ANOVA no existe interacción.
modelo1b <- aov(D1$Biomasa1~D1$Dosis_insecticida1+D1$Numero_aplicaciones1)
TukeyHSD(modelo1b)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = D1$Biomasa1 ~ D1$Dosis_insecticida1 + D1$Numero_aplicaciones1)
##
## $`D1$Dosis_insecticida1`
## diff lwr upr p adj
## 3.75-3.25 0.27263391 -0.17270866 0.7179765 0.2738312
## 4.25-3.25 0.35602957 -0.08931299 0.8013721 0.1261203
## 4.25-3.75 0.08339567 -0.36194690 0.5287382 0.8751713
##
## $`D1$Numero_aplicaciones1`
## diff lwr upr p adj
## 5-4 -0.24371366 -0.6890562 0.2016289 0.3479341
## 6-4 -0.32456967 -0.7699122 0.1207729 0.1711335
## 6-5 -0.08085601 -0.5261986 0.3644866 0.8821270Comparando lasmedias del efecto de las dosis del insecticida sobre la biomasa del cultivo se observa que el p-valor ajustado en todos los casos es mayor a 0.05, por lo que usar una doses o otra es lo mismo, se supone tendría el mismo efecto, es decir, todas las dosis son iguales entre sí ya que existe evidencia estadística para decir que tienen el mismo efecto sobre la biomasa. Aquí ya vamos contradiciendo el análisis descriptivo hecho anteriormente, en el que con la dosis 4.25 las plantas presentan mayor biomasa. Por otro lado, en las comparaciones de medias del efecto del número de aplicaciones sobre la biomasa, estas tienen el mismo comportamiento que las dosis del insecticida y por tanto, su efecto sobre la biomasa es igual.
Revisión de supuestos
- Prueba de normalidad de los residuales \[H_O= Residuales~normales\] \[H_O= Residuales~no~normales\]
# Prueba de normalidad de los residuales
residuales=modelo1b$residuals; residuales## 1 2 3 4 5 6
## -0.48411735 -0.14761743 -0.45238058 -0.11607671 0.03207085 -0.25891959
## 7 8 9 10 11 12
## -0.09802343 -0.03590285 -0.05094965 0.41038930 0.10233013 0.23758720
## 13 14 15 16 17 18
## 0.21133006 0.25860398 -0.01651421 0.44641666 -0.07578647 0.03756007shapiro.test(residuales)##
## Shapiro-Wilk normality test
##
## data: residuales
## W = 0.96635, p-value = 0.727Como el p-valor 0.06696>0.05, los residuales son normales por tanto se cumple el supuesto de normalidad.
- Prueba de homocedasticidad de residuales
\[H_0= Varianzas~iguales\] \[H_0= Al menos~una~varianza~diferente\]
# Prueba de homocedasticidad de residuales dosis
HR=bartlett.test(residuales~D1$Dosis_insecticida1)
HR##
## Bartlett test of homogeneity of variances
##
## data: residuales by D1$Dosis_insecticida1
## Bartlett's K-squared = 2.2328, df = 2, p-value = 0.3275El p-valor 0.2726>0.05 por tanto, se cumple la igualdad de varianzas o supuesto de homocedasticidad.
HR2=bartlett.test(residuales~D1$Numero_aplicaciones1)
HR2##
## Bartlett test of homogeneity of variances
##
## data: residuales by D1$Numero_aplicaciones1
## Bartlett's K-squared = 0.39774, df = 2, p-value = 0.8197El p-valor 0.7646>0.05 por tanto, se cumple la igualdad de varianzas o supuesto de homocedasticidad.
- Supuesto de independencia de residuales
# Grafico de independencia de residuales
plot(residuales, pch=16, col="red", main="Gráfico de residuales")
En el gráfico no se observa tendencia.
Ahora visualizaremos la tabla de medias, para comparar nuevamente el efecto de los tratamientos sobre el crecimiento del cultivo (biomasa):
t <- tapply(Biomasa1, list(Dosis_insecticida1, Numero_aplicaciones1), mean)
addmargins(t, FUN= mean)## 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.940436De la tabla anteior, se puede concluir que el mejor tratamiento, o mejor dicho, el que afecta menos el crecimiento del cultivo (mayor biomasa) si es la dosis 4,25 con 4 aplicaciones de acuerdo a las margenes y sin margenes también, estos quiere decir que no hay interacción, exactamente lo mismo que habiamos concluido del análisis descriptivo.
- Gráfico de interacción
library(dplyr)## Warning: package 'dplyr' was built under R version 4.0.3##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2)
D1 %>%
group_by(Dosis_insecticida1,Numero_aplicaciones1) %>%
summarise(medias_biomasa = mean(Biomasa1)) -> mediasDQ## `summarise()` regrouping output by 'Dosis_insecticida1' (override with `.groups` argument)ggplot(mediasDQ)+
aes(x=Dosis_insecticida1,y=medias_biomasa,col = Numero_aplicaciones1) +
geom_line(aes(group = Numero_aplicaciones1))+
geom_point()
A pesar de que en el ANOVA y en la comparación de medias habiamos concluido que no había interacción entre los factores Dosis del insecticida y Número de aplicaciones, en este gráfico se observa lo contrario, según este entre las aplicaciones 4 y 5 interaccionan (no en gran medida) con las 6 aplicaciones, pero por recomendación de profesor a veces los gráficos pueden ser engañosos y debemos confiar en los resultados del ANOVA y comparación de medias, en coclusión no hay interacción entre los dos 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?
- Covariable
set.seed(2020)
arcilla<-sort.default(runif(18, 0.2, 0.4), decreasing = TRUE) #covarable
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.2005165ARC=data.frame(D, arcilla)
ARC## Dosis_insecticida Numero_aplicaciones Repeticiones biomasa arcilla
## 2 3.75 4 r1 2.708826 0.3921445
## 10 3.25 4 r2 2.772692 0.3652331
## 16 3.25 6 r2 2.143359 0.3528828
## 4 3.25 5 r1 2.560519 0.3487672
## 13 3.25 5 r2 2.708666 0.3307115
## 6 4.25 5 r1 2.773705 0.3293806
## 17 3.75 6 r2 2.770350 0.3240412
## 8 3.75 6 r1 2.832470 0.3237004
## 14 3.75 5 r2 2.898280 0.3079385
## 5 3.75 5 r1 3.359619 0.2953782
## 9 4.25 6 r1 3.054099 0.2845458
## 1 3.25 4 r1 3.157896 0.2818575
## 12 4.25 4 r2 3.487669 0.2788452
## 11 3.75 4 r2 3.451547 0.2786236
## 15 4.25 5 r2 3.016111 0.2272194
## 7 3.25 6 r1 3.042156 0.2258305
## 3 4.25 4 r1 3.200552 0.2134769
## 18 4.25 6 r2 2.989329 0.2005165- Modelo del diseño incluyendo la covariable:
\[y_{ijk}=\mu+\alpha_i+\beta_j+(\alpha\beta)_{ij}+\delta(x_{ijk}-\bar{x}+\epsilon_{ijk}\\i:1\cdots3\\j:1,3\\k:1\cdots6\] Donde: - - \(y_{ijk}\) Es la \(ikj-esima\) observación de biomasa para el \(i-esimo\) nivel del factor 1 (Dosis_insecticida) y el \(j-esimo\) nivel del factor 2 (Numero_aplicaciones) y la \(k-esima\) repetición.
- \(\mu\) Es la media general de Biomasa
- \(\alpha_i\) Es el efecto del \(i-esimo\) nivel del factor 1 (Dosis_insecticida).
- \(\beta_j\) Es el efecto del \(j-esimo\) nivel del factor 2 (Numero_aplicaciones).
- \((\alpha\beta)_{ij}\) Es la interacción del \(i-esimo\) nivel del factor 1 (Dosis_insecticida) con el \(j-esimo\) nivel del factor 2 (Numero_aplicaciones). -\(x_{ijk}\) es la medida de la covariable que se hace para \(y_{ijk}\)
- \(\bar{x}\) es la media de los valores de \(x_{ijk}\)
- \(\delta\) es el coeficiente de regresión que relaciona \(y_{ijk}\) con la covariable \(x_{ijk}\)
- \(\epsilon_{ijk}\) Es el error residual para el \(i-esimo\) nivel del factor 1 (Dosis_insecticida), \(j-esimo\) nivel del factor 2 (Numero_aplicaciones) y la \(k-esima\) repetición.
Ancova
- Se corre el modelo con la interacción
ancova1<-aov(ARC$biomasa~ARC$arcilla+(ARC$Dosis_insecticida*ARC$Numero_aplicaciones))
summary(ancova1)## Df Sum Sq Mean Sq F value Pr(>F)
## ARC$arcilla 1 0.6594 0.6594 12.238 0.00393
## ARC$Dosis_insecticida 1 0.0681 0.0681 1.264 0.28119
## ARC$Numero_aplicaciones 1 0.4277 0.4277 7.936 0.01454
## ARC$Dosis_insecticida:ARC$Numero_aplicaciones 1 0.0123 0.0123 0.229 0.64044
## Residuals 13 0.7005 0.0539
##
## ARC$arcilla **
## ARC$Dosis_insecticida
## ARC$Numero_aplicaciones *
## ARC$Dosis_insecticida:ARC$Numero_aplicaciones
## Residuals
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1En el análisis de covarianza se observa que la arcilla tiene un efecto significativo sobre la biomasa, el número de aplicaciones también presenta efecto sobre esta, esto se detecto ya que al incluir la covariable e el análisis disminuyen los grados de libertad. Es decir, se mejoró la exactitud del modelo reduciendo el error.
Adicionalmente, en este análisis tampoco hubo interacción de los factores (dosis y aplicaciones).
- Ahora se corre sin la interacción (lo mejor)
ancova2<-aov(ARC$biomasa~ARC$arcilla+ARC$Dosis_insecticida+ARC$Numero_aplicaciones)
summary(ancova2)## Df Sum Sq Mean Sq F value Pr(>F)
## ARC$arcilla 1 0.6594 0.6594 12.951 0.00291 **
## ARC$Dosis_insecticida 1 0.0681 0.0681 1.338 0.26676
## ARC$Numero_aplicaciones 1 0.4277 0.4277 8.399 0.01168 *
## Residuals 14 0.7128 0.0509
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Aquí el p.valor de la arcilla 0.00291<0.05, lo que indica que la cantidad de arcilla si esta relacionadasignificativamente con la variable respuesta “Biomasa”. También se encontró significancia en el Número de aplicaciones sobre la biomasa.
Supuestos del ancova - Supuesto de linealidad entre la covariable y la biomasa.
plot(ARC$biomasa~ARC$arcilla)
Q <-lm(ARC$biomasa~ARC$arcilla,data = D)
abline(Q, col="red")
Se observa que los datos se distribuyen cerca de la linea, cumplen el supuesto.
- Homogeneidad de las pendientes de regresión o interacción nula entre la covariable y la biomasa.
library(rstatix)
ARC %>% anova_test(biomasa ~ arcilla *(Dosis_insecticida+Numero_aplicaciones)+arcilla*Dosis_insecticida*Numero_aplicaciones)## Coefficient covariances computed by hccm()## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05
## 1 arcilla 1 10 11.401 0.007 *
## 2 Dosis_insecticida 1 10 1.146 0.310
## 3 Numero_aplicaciones 1 10 12.090 0.006 *
## 4 arcilla:Dosis_insecticida 1 10 6.960 0.025 *
## 5 arcilla:Numero_aplicaciones 1 10 0.204 0.661
## 6 Dosis_insecticida:Numero_aplicaciones 1 10 1.135 0.312
## 7 arcilla:Dosis_insecticida:Numero_aplicaciones 1 10 0.032 0.862
## ges
## 1 0.533
## 2 0.103
## 3 0.547
## 4 0.410
## 5 0.020
## 6 0.102
## 7 0.003De acuerdo a la tabla, se observa no hay interacción entre la covariable y el numero de aplicaciones pero si con la dosis ya que presenta un p-valor 0.025<0.05. No se cumple este supuesto.
- Normalidad de los residuos.
G<- lm(ARC$biomasa~ARC$arcilla)
summary(G)##
## Call:
## lm(formula = ARC$biomasa ~ ARC$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 ***
## ARC$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.00932shapiro.test(G$residuals)##
## Shapiro-Wilk normality test
##
## data: G$residuals
## W = 0.94972, p-value = 0.4206Se acepta el supuesto porqie el p.valor 0.4206>0.05.
- Homogeneidad de la varianza.
ARC$Trt <- interaction(ARC$Dosis_insecticida, ARC$Numero_aplicaciones)
bartlett.test(ancova2$residuals~ARC$Trt)##
## Bartlett test of homogeneity of variances
##
## data: ancova2$residuals by ARC$Trt
## Bartlett's K-squared = 5.8077, df = 8, p-value = 0.6688Se cumple el supuesto puesto que el p-valor 0.368>0.05.
- Independencia de los residuales
plot(G$residuals)
Se cumple este supuesto.
Gráfico de interacción
Dosis_insecticida2= as.factor(ARC$Dosis_insecticida)
Aplicaciones2 = as.factor(ARC$Numero_aplicaciones)
Biomasa2 = as.vector(ARC$biomasa)
D_2 <- data.frame(Dosis_insecticida2, Aplicaciones2, Biomasa2)
str(D_2)## 'data.frame': 18 obs. of 3 variables:
## $ Dosis_insecticida2: Factor w/ 3 levels "3.25","3.75",..: 2 1 1 1 1 3 2 2 2 2 ...
## $ Aplicaciones2 : Factor w/ 3 levels "4","5","6": 1 1 3 2 2 2 3 3 2 2 ...
## $ Biomasa2 : num 2.71 2.77 2.14 2.56 2.71 ...library(dplyr)
library(ggplot2)
D_2 %>%
group_by(Dosis_insecticida2, Aplicaciones2) %>%
summarise(medias_groups = mean(Biomasa2)) -> MEDIAS_MD## `summarise()` regrouping output by 'Dosis_insecticida2' (override with `.groups` argument)ggplot(MEDIAS_MD)+
aes(x=Dosis_insecticida2,y=medias_groups,col = Aplicaciones2) +
geom_line(aes(group = Aplicaciones2))+
geom_point()
Se obtiene un gráfico de interacción muy similar al anterior.
PUNTO 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.
Tabla ANOVA
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.
- Los datos se muestran a continuación:
library(readxl)
papa <- read_excel("C:/Users/CK0001/Desktop/ING.AGRONOMICA.UNAL/semestre 6/D.E/Parcial/papa.xlsx")
P = data.frame(papa) #Crear data.frame
P## parcela Finca variedad test rendimiento
## 1 1 1 1 1 9.76
## 2 1 1 1 2 9.24
## 3 1 1 2 1 11.91
## 4 1 2 1 1 9.02
## 5 2 1 1 1 10.65
## 6 2 1 1 2 7.77
## 7 2 1 2 1 10.00
## 8 2 2 1 1 13.69
## 9 3 1 1 1 6.50
## 10 3 1 1 2 6.26
## 11 3 1 2 1 8.02
## 12 3 2 1 1 7.95
## 13 4 1 1 1 8.08
## 14 4 1 1 2 5.28
## 15 4 1 2 1 9.15
## 16 4 2 1 1 7.46
## 17 5 1 1 1 7.84
## 18 5 1 1 2 5.91
## 19 5 1 2 1 7.43
## 20 5 2 1 1 6.11
## 21 6 1 1 1 9.00
## 22 6 1 1 2 8.38
## 23 6 1 2 1 7.01
## 24 6 2 1 1 8.58
## 25 7 1 1 1 12.81
## 26 7 1 1 2 13.58
## 27 7 1 2 1 11.13
## 28 7 2 1 1 10.00
## 29 8 1 1 1 10.62
## 30 8 1 1 2 11.71
## 31 8 1 2 1 14.07
## 32 8 2 1 1 14.56
## 33 9 1 1 1 4.88
## 34 9 1 1 2 4.96
## 35 9 1 2 1 4.08
## 36 9 2 1 1 4.76
## 37 10 1 1 1 9.38
## 38 10 1 1 2 8.02
## 39 10 1 2 1 6.73
## 40 10 2 1 1 6.99
## 41 11 1 1 1 5.91
## 42 11 1 1 2 5.79
## 43 11 1 2 1 6.59
## 44 11 2 1 1 6.55
## 45 12 1 1 1 7.19
## 46 12 1 1 2 7.22
## 47 12 1 2 1 5.77
## 48 12 2 1 1 8.33
## 49 13 1 1 1 7.93
## 50 13 1 1 2 6.48
## 51 13 1 2 1 8.12
## 52 13 2 1 1 7.43
## 53 14 1 1 1 3.70
## 54 14 1 1 2 2.86
## 55 14 1 2 1 3.95
## 56 14 2 1 1 5.92
## 57 15 1 1 1 4.64
## 58 15 1 1 2 5.70
## 59 15 1 2 1 5.96
## 60 15 2 1 1 5.88
## 61 16 1 1 1 5.94
## 62 16 1 1 2 6.28
## 63 16 1 2 1 4.18
## 64 16 2 1 1 5.24
## 65 17 1 1 1 9.50
## 66 17 1 1 2 8.00
## 67 17 1 2 1 11.25
## 68 17 2 1 1 11.14
## 69 18 1 1 1 10.93
## 70 18 1 1 2 12.16
## 71 18 1 2 1 9.51
## 72 18 2 1 1 12.71
## 73 19 1 1 1 11.95
## 74 19 1 1 2 10.58
## 75 19 1 2 1 16.79
## 76 19 2 1 1 13.08
## 77 20 1 1 1 4.34
## 78 20 1 1 2 5.45
## 79 20 1 2 1 7.51
## 80 20 2 1 1 5.21- Arbol que representa el diseño anidado ecalonado:
library(collapsibleTree)
collapsibleTree(P,hierarchy=c("Finca","variedad","test"))Para estimar la varianza se va a usar el método REML, para ello se usa la función lmer del paquete lm4:
library(lme4)## Warning: package 'lme4' was built under R version 4.0.3## 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 carmodelo_papa <- lmer(rendimiento ~ 1 + (1|parcela) + (1|parcela:Finca)+
(1|parcela:Finca:variedad), data = P)## boundary (singular) fit: see ?isSingularsummary(modelo_papa)## Linear mixed model fit by REML ['lmerMod']
## Formula: rendimiento ~ 1 + (1 | parcela) + (1 | parcela:Finca) + (1 |
## parcela:Finca:variedad)
## Data: P
##
## 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.
## parcela:Finca:variedad (Intercept) 1.2305 1.1093
## parcela:Finca (Intercept) 0.0000 0.0000
## parcela (Intercept) 7.0127 2.6481
## Residual 0.8795 0.9378
## Number of obs: 80, groups:
## parcela:Finca:variedad, 60; parcela:Finca, 40; parcela, 20
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 8.2369 0.6188 13.31
## convergence code: 0
## boundary (singular) fit: see ?isSingularSe observa que la varianza para la finca dentro de la parcela es 0, por lo tanto, no hay varianza atribuible a las parcelas anidadas en las fincas.
Porcentaje de la variación de las Parcelas y de las variedades dentro de la finca:
variacion_Parcelas= (100*7.0186)/(7.0186+1.2361+0.8743)
variacion_Parcelas## [1] 76.88246variacion_V= (100*1.2361)/(1.2361+7.0186+0.8743)
variacion_V## [1] 13.54037De acuerdos a estos resultados, se concluye que el 76.9% de la variación total se debe a la variabilidad entre parcelas y el 13,54% se debe a la variabilidad de las variedades dentro de la finca y la variabilidad dentro de parcela es insignificante (0). Esa gran variabilidad de las parcelas, agronómicamente puede significar que si los agricultores quieren obtener un rendimiento de papa menos variable, en primer lugar deben reducir la variación de las parcelas, esta puede darse porlas propiedades físico-quimicas y biológicas del suelo.
PUNTO 6
se tienen unos datos de potasio de muestras de suelos medidas en 8 diferentes laboratorios. Los datos son:
library(asbio)## Warning: package 'asbio' was built under R version 4.0.3## Loading required package: tcltkdata(K)
K_suelo<-K
K_suelo## K lab
## 1 296 B
## 2 260 B
## 3 341 B
## 4 359 B
## 5 323 B
## 6 321 B
## 7 287 B
## 8 413 B
## 9 335 B
## 10 315 D
## 11 330 D
## 12 326 D
## 13 354 D
## 14 266 D
## 15 348 D
## 16 343 D
## 17 284 D
## 18 324 D
## 19 351 E
## 20 302 E
## 21 395 E
## 22 357 E
## 23 400 E
## 24 187 E
## 25 376 E
## 26 283 E
## 27 198 E
## 28 327 F
## 29 354 F
## 30 308 F
## 31 274 F
## 32 324 F
## 33 305 F
## 34 347 F
## 35 297 F
## 36 305 F
## 37 326 G
## 38 301 G
## 39 316 G
## 40 312 G
## 41 297 G
## 42 280 G
## 43 300 G
## 44 319 G
## 45 286 G
## 46 218 H
## 47 280 H
## 48 241 H
## 49 226 H
## 50 243 H
## 51 199 H
## 52 205 H
## 53 225 H
## 54 227 H
## 55 338 I
## 56 303 I
## 57 341 I
## 58 311 I
## 59 355 I
## 60 269 I
## 61 284 I
## 62 279 I
## 63 339 I
## 64 359 J
## 65 318 J
## 66 313 J
## 67 352 J
## 68 334 J
## 69 356 J
## 70 342 J
## 71 299 J
## 72 353 JCompare 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.
- En primer lugar se realizará un diagrama de cajas para visualizar el comportamiento de nuestros datos y tener una idea inicial de cual laboratorio cuenta con datos menos variables o más confiables:
library(ggplot2)
ggplot(K_suelo, aes(x = lab, y = K))+
geom_boxplot()
Entonces, de acuerdo al diagrama de cajas, desde mi perspectiva observo que el laboratorio H pese a que presenta los niveles de K registrados más bajitos, estos no son tan variables, sin embargo por ser el más bajito teniendo en cuenta que las muestras enviadas a todos los laboratorios son de la misma finca, este no es confiable ya que los demás laboratorios presentan valores cercanos entre sí y solo el H es la excepción, y adicionalmente se resalta que presenta un dato atípico.
En cuanto al laboratorio que registran los niveles más altos de K, este es el E, pero a su vez es el laboratorio con más variabilidad entre sus datos, cosa que tampoco da buena referencia. Por ahora, de acuerdo a este diagrama el laboratorio J, para mí o desde mi punto de vista, es el que presenta datos altos y confiables de los niveles de K en el suelo de la finca, dado que sus datos están entre los más altos, no hay atípicos y presenta poca variablidad.
El diagrama de violin resulta muy útil para observar la forma de la distribución o como se concentran los niveles de K registrados por cada laboratorio:
ggplot(data = K_suelo , aes(x = lab , y = K)) +
geom_jitter(size = 1, color = 'gray', alpha = 0.5) +
geom_violin(aes(fill = lab), color = 'black', alpha = .8 ) +
geom_boxplot(color = 'black', alpha = 0.7) +
stat_summary(fun=median, geom="point", shape=23, size=2, color = "red")+
xlab('Laboratorio') +
ylab('K') +
ggtitle('Registro de K en los laboratorios') +
theme_minimal()
De este gráfico se sigue observando que el laboratorio J presenta altos niveles de K registrados pero el laboratorio que tiene datos más densos o cncentrados es el G y si nos devolvemos al diagrama de cajas efectivamente el G es el que presenta menos variabilidad de sus datos. También se comprueba que el laboratorio E es el más variable (menos denso). Hata aquí se pueden consigerar los laboratorios G y J como los mejores, sin embargo para tomar una desición de cual prueba recomendaría, este análisis descriptivo lo debemos complementar con análisis de medidas y tablas.
-Primero haré una comparación mediante la tabla de medias:
medias_K <- tapply(K_suelo$K, list(K_suelo$lab), FUN = mean)
medias_K## B D E F G H I J
## 326.1111 321.1111 316.5556 315.6667 304.1111 229.3333 313.2222 336.2222- Lo anteior complementado con el coeficiente de variación teniendo en cuenta que un CV mayor >20% no es recomendable, nos ayuda a definir el mejor laboratorio, entonces:
desviacion=tapply(K_suelo$K, list(K_suelo$lab), sd)
coef_var = 100 * desviacion/medias_K
coef_var## B D E F G H I J
## 13.616722 9.110586 25.382896 7.935570 5.055402 10.417315 10.023481 6.426969De acuerdo al diagrama de violin, tenia como pocibles opciones de mejor laboratorios el G y J, ya con a la tabla de medias puedo ver que el J sobresale un poco (mayor media) y finalmente el CV me dice que el laboratorio G es el que presenta menor variación pero la diferencia con el J no es mucha, eso quiere decir que si me toca recomendar un laboratorio para el analisis de K en el suelo, recomendaría en laboratorio J en primer lugar y el G en egundo lugar.
Sabiendo que son muestras mezcladas de una misma finca, ¿ Se perciben diferencias en las medias 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).
- Enfoque no paramétrico: Prueba de Krustall-Wallis: Este es un test que emplea los rangos para contrastar la hipótesis de que k muestras han sido obtenidas de una misma población. Entonces se realizará un análisis de varianza de los rangos,compararando los promedios de los estos para cada laboratorio.
Hipótesis: \[H_0= \mu_{rango}~es~igual~para~todos~los~laboratorios \] \[H_0= \mu_{rango}~es~diferente~en~al~menos~un~laboratorio\]
- Se calcula la media de los rangos de los datos de los laboratorios y se realizará un análisis gráfico:
which.min(K_suelo$K)## [1] 24K_rangos <-rank(K_suelo$K)
tapply(K_rangos, K_suelo$lab, mean)## B D E F G H I J
## 42.444444 42.444444 42.777778 38.277778 30.611111 7.611111 36.888889 50.944444Con esta tabla ya se puede ir concluyendo que los rangos difieren, principalmente el del laboratorio H.
boxplot(K_rangos~K_suelo$lab)
Con el diagrama de cajas se observa que todos los laboratorios excepto el H mantienen un comportamiento similar. El laboratorio H presenta diferencia significativas.
Ahora, se realiza la prueba Krustall-Wallis:
kruskal.test(K_suelo$K~K_suelo$lab)##
## Kruskal-Wallis rank sum test
##
## data: K_suelo$K by K_suelo$lab
## Kruskal-Wallis chi-squared = 24.482, df = 7, p-value = 0.000937Evidentemente se rechaza la hipotesis nula (p-valor 0.000937<0.05), por lo que acepta la hipotesis alterna que dice la nedia de los rangos es diferente en al menos un laboratorio y de acuerdo al análisis gráfico, es diferente significativamente en el laboratorio H.
Finalmente, con la prueba de Nemenyi se pueden comparar las medias por pares de laboratorio y encontrar las diferencias.
PMCMR::posthoc.kruskal.nemenyi.test(K_suelo$K~K_suelo$lab)## 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_suelo$K by K_suelo$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: noneCon esta prueba se observa que en la mayoría de parejas se obtiene un p-valor mayor a 0.05. Sin embargo, por otra parte, los laboratorios presentan diferenciassignificativas estadisticamente pero únicamente con el laboratorio H , pues al compararse los otros, este presentó un p-valor menor de 0.05 en la mayoría. Por lo tanto, definitivamente se rechaza la hipótesis nula que indica que todos los laboratorios son iguales.
PUNTO 7. Strip plot design
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 en franjas
Se va a analizar el rendimiento de frijol utilizando 3 variedades con 3 niveles de aplicaciones de urea durante el ciclo productivo. Este experimento se dispone en 6 bloques de 3 parcelas principales (variedad) y estas a su vez se dividen en 3 subparcelas (número de aplicaiones de urea). Se busca encontrar con cual tratamiento principalmente con cual dosis se obtiene un mayor rendimiento. Entonces los datos fueron tomados de [https://sscnars.icar.gov.in/R_manual/R-Scripts%20for%20Statistical%20analysis%20using%20R-studio.pdf] y se creó un excel con ellos, estos son:
library(readxl)
Frijol <- read_excel("C:/Users/CK0001/Desktop/ING.AGRONOMICA.UNAL/semestre 6/D.E/Parcial/frijol.xlsx")
Frijol## # A tibble: 54 x 4
## Bloque variedad_frijol aplicaciones_Urea rendimiento_frijol
## <dbl> <dbl> <dbl> <dbl>
## 1 1 1 1 104.
## 2 1 1 2 126.
## 3 1 1 3 115.
## 4 1 2 1 104.
## 5 1 2 2 119.
## 6 1 2 3 106.
## 7 1 3 1 138.
## 8 1 3 2 133.
## 9 1 3 3 129.
## 10 2 1 1 91.6
## # ... with 44 more rowsFrijol_DQ = as.data.frame(Frijol)
Frijol_DQ## Bloque variedad_frijol aplicaciones_Urea rendimiento_frijol
## 1 1 1 1 103.60
## 2 1 1 2 126.25
## 3 1 1 3 115.39
## 4 1 2 1 103.80
## 5 1 2 2 119.32
## 6 1 2 3 106.40
## 7 1 3 1 138.30
## 8 1 3 2 132.81
## 9 1 3 3 128.90
## 10 2 1 1 91.61
## 11 2 1 2 103.30
## 12 2 1 3 101.56
## 13 2 2 1 116.65
## 14 2 2 2 125.40
## 15 2 2 3 144.65
## 16 2 3 1 106.60
## 17 2 3 2 129.35
## 18 2 3 3 138.30
## 19 3 1 1 126.91
## 20 3 1 2 119.90
## 21 3 1 3 134.19
## 22 3 2 1 131.19
## 23 3 2 2 119.73
## 24 3 2 3 127.11
## 25 3 3 1 125.80
## 26 3 3 2 113.80
## 27 3 3 3 125.35
## 28 4 1 1 129.61
## 29 4 1 2 118.39
## 30 4 1 3 123.00
## 31 4 2 1 120.81
## 32 4 2 2 120.26
## 33 4 2 3 125.60
## 34 4 3 1 104.71
## 35 4 3 2 135.60
## 36 4 3 3 114.31
## 37 5 1 1 128.29
## 38 5 1 2 129.40
## 39 5 1 3 123.71
## 40 5 2 1 131.29
## 41 5 2 2 122.85
## 42 5 2 3 138.20
## 43 5 3 1 123.31
## 44 5 3 2 142.71
## 45 5 3 3 152.21
## 46 6 1 1 121.06
## 47 6 1 2 112.95
## 48 6 1 3 115.91
## 49 6 2 1 119.60
## 50 6 2 2 122.45
## 51 6 2 3 126.61
## 52 6 3 1 126.96
## 53 6 3 2 130.94
## 54 6 3 3 131.70Modelo del diseño
\[y_{ijk}=\mu+\alpha_i+\gamma_k+\eta_{ik}+\beta_j+\ (\alpha\beta)_{ij} +\epsilon_{ijk\\i:1\cdots 3\\j:1,3\\k:1\cdots 6}\\Condiciones~Laterales~Respectivas\]
Donde:
\(Y_{ijk}\) es el rendimiento del frijol
\(\mu\) es la media general
\(\alpha_i\) es el efecto fijo de la variedad
\(\gamma_{ik}\) es el efecto fijo del bloqueo
\(\eta_{ik}\) es el error en parcela principal
\(\beta_j\) es el efecto fijo de numero de aplicaiones de urea
\((\alpha\beta)_{ij}\) es el efecto de la interacción entre ambos factores
\(\epsilon_{ijk}\) es el error en la subparcela
str(Frijol_DQ)## 'data.frame': 54 obs. of 4 variables:
## $ Bloque : num 1 1 1 1 1 1 1 1 1 2 ...
## $ variedad_frijol : num 1 1 1 2 2 2 3 3 3 1 ...
## $ aplicaciones_Urea : num 1 2 3 1 2 3 1 2 3 1 ...
## $ rendimiento_frijol: num 104 126 115 104 119 ...variedad_frijol2 = factor(Frijol_DQ$variedad_frijol)
aplicacionesU = factor(Frijol_DQ$aplicaciones_Urea)
Bloqueo2 = as.numeric(Frijol_DQ$Bloque)
Rendimiento_frijol=as.vector(Frijol_DQ$rendimiento_frijol)
Frijol_2020=data.frame(variedad_frijol2,aplicacionesU,Rendimiento_frijol)
str(Frijol_2020)## 'data.frame': 54 obs. of 3 variables:
## $ variedad_frijol2 : Factor w/ 3 levels "1","2","3": 1 1 1 2 2 2 3 3 3 1 ...
## $ aplicacionesU : Factor w/ 3 levels "1","2","3": 1 2 3 1 2 3 1 2 3 1 ...
## $ Rendimiento_frijol: num 104 126 115 104 119 ...Ahora un analisis de varianza ANOVA.
attach(Frijol_2020)## The following objects are masked _by_ .GlobalEnv:
##
## aplicacionesU, Rendimiento_frijol, variedad_frijol2library(agricolae)## Warning: package 'agricolae' was built under R version 4.0.3modelo_Frijol = with(Frijol_2020, strip.plot(Bloqueo2, variedad_frijol2, aplicacionesU, Rendimiento_frijol))##
## ANALYSIS STRIP PLOT: Rendimiento_frijol
## Class level information
##
## variedad_frijol2 : 1 2 3
## aplicacionesU : 1 2 3
## Bloqueo2 : 1 2 3 4 5 6
##
## Number of observations: 54
##
## model Y: Rendimiento_frijol ~ Bloqueo2 + variedad_frijol2 + Ea + aplicacionesU + Eb + aplicacionesU:variedad_frijol2 + Ec
##
## Analysis of Variance Table
##
## Response: Rendimiento_frijol
## Df Sum Sq Mean Sq F value Pr(>F)
## Bloqueo2 5 1246.19 249.24 4.0447 0.010651 *
## variedad_frijol2 2 869.34 434.67 1.9445 0.193481
## Ea 10 2235.36 223.54 3.6276 0.006836 **
## aplicacionesU 2 427.31 213.66 1.9989 0.186087
## Eb 10 1068.90 106.89 1.7346 0.141300
## aplicacionesU:variedad_frijol2 4 219.64 54.91 0.8911 0.487390
## Ec 20 1232.42 61.62
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## cv(a) = 12.1 %, cv(b) = 8.4 %, cv(c) = 6.4 %, Mean = 123.1224De acuerdo al ANOVA obtenido, se observa que que no existe interacción entre ambos factores (variedades y número de aplicaciones de urea). Por otra parte, ninguno de los dos factores tienen un efecto significativo en el rendimiento del frijol. Sin embargo, también se obrserva que los bloques con su p-valor 0.010651<0.05 difieren significativamente, esto no debería ser así.
- Análisis descriptivo para establecer el posible mejor tratamiento:
lattice::bwplot(Frijol_2020$Rendimiento_frijol~Frijol_2020$aplicacionesU|Frijol_2020$variedad_frijol2)
De acuerdo al gráfico, se alcanza un mayor rendimiento del frijol con dos aplicaciones de urea con la variedad 3.
Para analizar mejor se visualizará la tabla de medias con margenes:
Q <- tapply(Frijol_2020$Rendimiento_frijol, list(Frijol_2020$variedad_frijol2, Frijol_2020$aplicacionesU), mean)
addmargins(Q, FUN= mean)## Margins computed over dimensions
## in the following order:
## 1:
## 2:## 1 2 3 mean
## 1 116.8467 118.3650 118.9600 118.0572
## 2 120.5567 121.6683 128.0950 123.4400
## 3 120.9467 130.8683 131.7950 127.8700
## mean 119.4500 123.6339 126.2833 123.1224De acuerdo a la tabla, sin las margenes el mejor rendimiento se da en la variedad 3 con 3 aplicaciones, y con las margenes, el mayor rendimiento se alcanza con la misma variedad 3 y las 3 aplicaciones, por lo que si quiero cuktuvar frijol, lo mejor es que siembre la variedad 3 y durante su ciclo productivo le haga 3 aplicaciones de urea.
- Añadir covariable: En este caso, se escogió como covariable el pH del suelo, debido a que este tiene gran influencia en el desarrollo y respectivo rendimiento de los cultivos.
set.seed(2020)
pH<-sort.int(runif(54, 4.7, 7.5), decreasing = TRUE) #covarable
pH## [1] 7.458451 7.399405 7.390023 7.376265 7.354340 7.338733 7.314803 7.227921
## [9] 7.177113 7.175829 7.119025 7.013264 6.991501 6.934219 6.840359 6.794666
## [17] 6.782740 6.661868 6.529961 6.511328 6.449555 6.436577 6.431805 6.320559
## [25] 6.300349 6.256451 6.230803 6.211139 6.165090 6.159959 6.035295 6.018364
## [33] 5.962661 5.883641 5.846005 5.803832 5.800730 5.728382 5.694284 5.507169
## [41] 5.444978 5.330924 5.321563 5.251028 5.164459 5.164089 5.081072 5.061627
## [49] 4.918039 4.900113 4.888676 4.730913 4.707232 4.703326Frijol_pH=data.frame(Frijol_DQ, pH)
Frijol_pH## Bloque variedad_frijol aplicaciones_Urea rendimiento_frijol pH
## 1 1 1 1 103.60 7.458451
## 2 1 1 2 126.25 7.399405
## 3 1 1 3 115.39 7.390023
## 4 1 2 1 103.80 7.376265
## 5 1 2 2 119.32 7.354340
## 6 1 2 3 106.40 7.338733
## 7 1 3 1 138.30 7.314803
## 8 1 3 2 132.81 7.227921
## 9 1 3 3 128.90 7.177113
## 10 2 1 1 91.61 7.175829
## 11 2 1 2 103.30 7.119025
## 12 2 1 3 101.56 7.013264
## 13 2 2 1 116.65 6.991501
## 14 2 2 2 125.40 6.934219
## 15 2 2 3 144.65 6.840359
## 16 2 3 1 106.60 6.794666
## 17 2 3 2 129.35 6.782740
## 18 2 3 3 138.30 6.661868
## 19 3 1 1 126.91 6.529961
## 20 3 1 2 119.90 6.511328
## 21 3 1 3 134.19 6.449555
## 22 3 2 1 131.19 6.436577
## 23 3 2 2 119.73 6.431805
## 24 3 2 3 127.11 6.320559
## 25 3 3 1 125.80 6.300349
## 26 3 3 2 113.80 6.256451
## 27 3 3 3 125.35 6.230803
## 28 4 1 1 129.61 6.211139
## 29 4 1 2 118.39 6.165090
## 30 4 1 3 123.00 6.159959
## 31 4 2 1 120.81 6.035295
## 32 4 2 2 120.26 6.018364
## 33 4 2 3 125.60 5.962661
## 34 4 3 1 104.71 5.883641
## 35 4 3 2 135.60 5.846005
## 36 4 3 3 114.31 5.803832
## 37 5 1 1 128.29 5.800730
## 38 5 1 2 129.40 5.728382
## 39 5 1 3 123.71 5.694284
## 40 5 2 1 131.29 5.507169
## 41 5 2 2 122.85 5.444978
## 42 5 2 3 138.20 5.330924
## 43 5 3 1 123.31 5.321563
## 44 5 3 2 142.71 5.251028
## 45 5 3 3 152.21 5.164459
## 46 6 1 1 121.06 5.164089
## 47 6 1 2 112.95 5.081072
## 48 6 1 3 115.91 5.061627
## 49 6 2 1 119.60 4.918039
## 50 6 2 2 122.45 4.900113
## 51 6 2 3 126.61 4.888676
## 52 6 3 1 126.96 4.730913
## 53 6 3 2 130.94 4.707232
## 54 6 3 3 131.70 4.703326ANCOVA -Primero se observa si vale la pena realizar un análisis ANCOVA para este modelo, entonces se analiza la covariable pH con el rendimiento del frijol, de manera gráfica:
plot(Frijol_pH$pH, Frijol_pH$rendimiento_frijol)
- Ahora se calcula el test de correlación.
cor.test(Frijol_pH$pH, Frijol_pH$rendimiento_frijol)##
## Pearson's product-moment correlation
##
## data: Frijol_pH$pH and Frijol_pH$rendimiento_frijol
## t = -2.451, df = 52, p-value = 0.01764
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.54280081 -0.05914897
## sample estimates:
## cor
## -0.3218127Al obtener un p_valor (0.01764) mayor al 5% se confirma que existe correlación entre ambas variables. Entonces se corre el ANCOVA:
ancova_frijol<-aov(Frijol_pH$rendimiento_frijol~Frijol_pH$pH+(Frijol_pH$variedad_frijol*Frijol_pH$aplicaciones_Urea))
summary(ancova_frijol)## Df Sum Sq Mean Sq F value
## Frijol_pH$pH 1 756 755.9 6.846
## Frijol_pH$variedad_frijol 1 642 642.2 5.817
## Frijol_pH$aplicaciones_Urea 1 377 376.8 3.413
## Frijol_pH$variedad_frijol:Frijol_pH$aplicaciones_Urea 1 114 113.9 1.032
## Residuals 49 5410 110.4
## Pr(>F)
## Frijol_pH$pH 0.0118 *
## Frijol_pH$variedad_frijol 0.0197 *
## Frijol_pH$aplicaciones_Urea 0.0707 .
## Frijol_pH$variedad_frijol:Frijol_pH$aplicaciones_Urea 0.3147
## Residuals
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1De acuerdo a los resultados, el pH al presentar un p-valor 0.0118 <0.05 influye significativamente en el rendimiento del frijol, al igual que la variedad. Se observa que no hay interacción entre los dos factores.
PUNTO 8
Realice un resumen con la nota que aparece en las siguientes direcciones sobre:
El uso de los diseños en parcelas divididas:
Inicialmente, cuando realizamos experimentos es muy importante comprender el diseño de este, ya que permite plantear estrategias para seleccionar, controlar, analizar e interpretar diferentes condiciones de un proceso de manera objetiva y sistemática, es decir, el diseño de experimentos permite mejorar o maximizar un proceso y la calidad de sus productos.
Entre estos diseños encontramos el experimento de parcelas divididas o Split plot Design, el cual es recomendable su uso cuando se tiene una situación donde los factores presentan niveles dificiles de cambiar o aleatorizar, estos factores que son difícil de aleatorizar completamente se debe a limitaciones de tiempo o costo. Entonces, en un experimento de parcelas divididas, los niveles del factor difícil de cambiar se mantienen constantes durante varias corridas experimentales, las cuales son tratadas como una parcela completa o parcela principal, mientras que los factores fáciles de cambiar varían en estas corridas, y cada combinación se considera una parcela subdividida dentro de la parcela completa es decir subparcelas dentro de la parcela principal.
Las parcelas divididas son originarias de entornos agrícolas por lo que se estas se refieren a parcelas de tierra donde por ejemplo la aplicación del mismo fertilizante se puede considerar la parcela principal y el tipo de semilla como las subparcelas porque estas se pueden aleatoriza de una manera más fácil. Este tipo de diseños también se puede dar en ámbitos industriales, y en este caso la parcela principal puede ser la temperatura de un horno, la configuración de equipos complejos e incluso la mezcla de ingredientes, estos factores difíciles de cambiar están ligados al presupuesto o al tiempo.
El número de parcelas completas corresponderá al número de veces que se repita la aplicación del factor de parcela completa, mientras que el número total de experimentos es equivalente al número de subparcelas.
Al analizar un diseño de parcelas divididas puede confundirlo con un diseño completamente al azar, pues puede darse un trazado dividido involuntario, donde se ejecuta un experimento de parcela dividida pero al analizarlo puede tratarse erróneamente de un diseño completamente aleatorio. Aquí todas las corridas se tratan como independientes lo cual pede llevar a conclusiones erróneas sobre los efectos de los factores, es decir, llegar a considerarlo significativo o no significativo cuando si es importante.
En cuanto al orden del experimento, al utlizar parcelas divididas el orden de las ejecuciones presenta una aleatorización más estructurada.Por otra parte los datos deben analizarse como un diseño de parcelas divididas lo cual nos llevará a conclusiones válidas con información precisa sobre los efectos. Entonces para el orden de ejecución de este tipo de diseño se realizan dos aleatorizaciones independientes: una para determinar el orden en el que se ejecutan las parcelas completas y determinar el orden en el que se ejecutan todas las gráficas; y recoger observaciones dentro de cada parcela otra para las observaciones dentro de cada parcela. A su vez, esto conducirá a dos términos de error en el modelo donde uno se atribuye a la parcela principal y el otro a la subparcela, y estos afectan a la forma como debe realizarse el análisis y como se comparan los diseños.
Si la estimación o la predicción son importantes, la parcela dividida involuntaria puede producir un diseño que es tan sólo 50% eficiente,es decir, este enfoque de “dejarlo al azar” tiene un alto potencial para decepcionar. Pero esta estimación del diseño completamente aleatorio puede superarse al elegir el experimento de parcela dividida correcto, este último puede producir estimaciones de parámetros más precisas.
En cuanto a costos, estos se pueden calcular como la suma ponderada del número de parcelas completas y el número total de subparcelas. Aquí se puede justificar este tipo de experimento ya que se puede seleccionar un diseño con menos parcelas enteras (ya que suelen ser costosos), lo que representaría un gran ahorro.
Al momento de legir un buen diseño, debemos tener en cuenta aspectos cualitativos como el equilibro entre el número de observaciones por parcela completa y el número de niveles de cada factor y cuantitativos como buena estimación de los términos de error y el rendimiento. Aquí es de gran importancia considerar el óptimo del experimento y los criterios para alcanzar el objetivo del proyecto. Los diseños de parcelas divididas, si se eligen bien pueden proporcionar mayor información al experimentador.
PUNTO 9
Seleccionar un artículo científico de una revista de agronomía donde se haya utilizado un diseño en parcelas divididas.
Revista:Agronomía Mesoamericana. 27(1):167-176. 2016
Artículo:Evauación de tres variedades de frijol Pinto bajo riego y sequía en Durango, México.
Autores:Aurelio Pedroza Sandoval, Ricardo Trejo Calzada, Ignacio Sánchez Cohen, José Alfredo Samaniego Gaxiola & Luis Gerardo Yáñez Chávez.
Hacerlas críticas constructivas sobre:
La mención de la estructura factorial
Se usó un diseño de bloques al azar en un arreglo de parcelas divididas con 3 repeticiones.
Las parcelas principales fueron los tratamientos de riego (contenido de humedad a capacidad de campo CC: 22 - 26%) y sequía (contenido de humedad a punto de marchitez permanente PMP: 16 – 20).
Las parcelas chicas o subparcelas fueron 3 variedades de frijol Pinto: Pinto Centauro, Pinto Americano y Pinto Saltillo.
De lo anterior, se tienen dos factores, el primer factor “A” presenta 2 niveles correspondientes a los tratamiento de riego y sequía. Por su parte el factor “B” tiene 3 niveles que son las 3 variedades de frijol Pinto. De este modo, puedo decir que se trata de un diseño con estructura factorial AxB, es decir, 2x3 (2 niveles de A por 3 niveles de B)para un total de 6 tratamientos con 3 replicaciones cada uno, por lo que son 18 observaciones en total.
Las variables morfométricas respuesta que analizaron fueron altura de planta, porcentaje de emergencia o germinación de plántulas respecto al total de semillas sembradas, porcentaje de plantas marchitas, cobertura vegetal y variables fisiológicas de fotosíntesis (F), conductancia, transpiración (Tr) y la eficiencia en el uso de agua (F/Tr).
La razón de colocar cada factor en la parcela principal o en la subparcela
Aunque ellos no especifican el porque escogieron las condiciones de riego y sequía como parcela principal y las variedades como subparcela, yo lo interpreto a partir de lo siguiente: “Cuando los niveles de un factor requieren de gran cantidad de material experimental frente a otros factores, como ejemplo podemos tener las siguientes situaciones: uso de riego, métodos de aplicación de fertilizantes, etc., serían más factibles usarlos como parcela principal que como subparcela y se elige como subparcelas los factores que se desea estudiar con mayor precisión” DISEÑO EN PARCELAS DIVIDIDAS.
De acuerdo a lo anterior, desde el objetivo que ellos se plantearon el cual fue evaluar el crecimiento y fisiología de tres variedades de frijol en condiciones de riego y sequía, vemos que el factor que ellos querían estudiar con mayor precisión eran las variedades, por eso las escogieron como subparcelas, además de que el factor riego y sequía exige parcelas relativamente grandes y por lo que sus niveles van a ser más difícil de aleatorizar.
La revisión de supuestos para el análisis de varianza
- No hicieron revisión de los supuestos del análisis de varianza en el artículo.
La tabla del análisis de varianza
¡No se observa tabla de ANOVA!.
Tal vez ellos la omitieron y mostraron directamente los resultados de la prueba tukey. Esto lo relaciono a que como no es suficiente con un análisis de varianza (ANOVA) para saber si en la comparativa de diferentes tratamientos (o experimentos) aplicada a varias muestras cumplen la hipótesis nula (Ho: “todos los tratamientos son iguales”) o por el contrario se cumple la hipótesis alternativa (Ha: “por lo menos uno de los tratamientos es diferente”), se aplica la prueba Tukey, como no muestran la tabla de analisis de varianza (me parece que es un gran error no incluirla) y aplican la prueba tukey se puede decir que el análisis de varianza que ellos hicieron rechazó la hipotesis nula, y es por ello que las únicas 3 tablas que se observan corresponden a las medias de las diferentes variables respuesta estudiadas para las tres variedades bajos las condiciones de riego y sequía, donde mediante la prueba Tukey (P≤0,05), las cifras con las mismas letras dentro de una misma columna, son estadísticamente iguales .
Algo importante es que el análisis estadístico fue en parcelas divididas en el tiempo, con el interés de conocer cuándo hubo o no un crecimiento significativo en el tiempo, independientemente de la condición de humedad en el suelo.
El uso de muchos análisis de varianzas en lugar de uno solo multivariante
- No mencionan cuantos ni muestran los resultados de los análisis de varianza, pero basada en las pruebas tukey que fueron una para cada tabla (porcentaje de emergencia de plántulas y de marchitez de plantas, Crecimiento promedio y Incremento de la cobertura vegetal promedio) asumo que son análisis de varianza para cada una de estas. Por otra parte, se observan dos conjuntos de gráficas que realmente me confunden porque parece que cada una tiene sus prueba tukey y por ende serián ANOVAs diferentes, pero el hecho de estar en conjunto no sé si signifique que ahi se realizó un anova multivariante.
El método de comparaciones de medias después del Anova
- Realizaron prueba de rango múltiple de medias Tukey, donde:
En la tabla 1, no hubo diferencia significativa en la nacencia de plántulas por efecto de variedad; en tanto que una vez emergida la plántula, se produjo una marchitez asociada a patógenos del suelo (Rhizoctonia solani y Fusarium solani), la cual fue significativamente mayor en la variedad Pinto Americano (16,0%) y moderadamente menor en las variedades Pinto Centauro y Pinto Saltillo (12,5% y 13,7%, respectivamente), sin diferencia estadística entre ambas.
La tabla dos y la tres son iguales, tal vez se equivocaron al colocarlas, sin embargo ellos alcanzan a mencionar que en las tres variedades de frijol evaluadas en este estudio, destacó la variedad Pinto Centauro, la cual siempre fue de mayor altura respecto a las otras dos variedades.
En la tabla 3, correspondiente a incremento de la cobertura vegetal promedio, destacó la variedad Pinto Centauro, un comportamiento intermedio fue observado en la variedad Pinto Americano y el de menor respuesta el Pinto Saltillo.
La interacción de factores
- De acuerdo a lo que mencionan los autores, exite interacción entre la variedad y el contenido de humedad, donde las tres variedades disminuyeron su crecimiento y desarrollo al pasar de la condición de riego a las de sequía, ya que a que el déficit hídrico disminuye significativamente el crecimiento, afectando de forma negativa la captación de radiación solar y por ende la producción de materia seca, debido a una menor conductancia y fotosíntesis.
La presencia de bloques
- Asumiendo que el número de repeticiones es igual al número de bloques, entonces hay 3 bloques al azar.
El balanceo o desbalanceo
- Asumo que es un diseño balanceado porque no se menciona pérdidas de datos.
La definición clara de la unidad experimental
- La siembra se estableció en surcos de 0,30m de altura y 0,8 m de ancho, con una distancia de 0,15 m entre plantas y 0,20m entre surcos. Se sembró a dos hileras en cada surco a una distancia de 0,30m y una cintilla de riego entre hileras. Cada unidad experimental fue de cuatro surcos de 6m de longitud cada uno, descartando los dos surcos extremos, los dos surcos medios correspondieron a la parcela útil, a partir de los cuales se seleccionaron cuatro plantas al azar y se midieron las respectivas variables.
Software utilizado y librería específica (en caso de ser R)
- Mediante uso del programa estadístico SAS, versión 9.0, se realizaron análisis estadísticos de varianza y prueba de rango múltiple de medias Tukey.
Otros aspectos que considere de interés
Este artíclo presenta falencias, principalmente por la falta de las tablas ANOVA, por otra parte, observando fijamente encontré que la tabla dos y la tres son exactamente iguales, los datos de ambas corresponden al incremento de la cobertura vegetal promedio \(cm^2\), cuando en realidad en la tabla dos deberían ir los datos del crecimiento promedio en \(cm\), es decir, que hubo una confusión con las gráficas y resultó faltando una “crecimiento”.