Parcial 1
Parcial 1 Diseño de Experimentos
Universidad Nacional de Colombia
Integrantes:
Erika Milena Aparicio Aponte
Ana María Montaño Hernández
1. Contenido de clorofila en 72 macetas.
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 883. Con esta información complete la tabla del ANOVA que se muestra a continuación:
Para la varianza se utiliza la cédula de Ana María Montaño con U=8
Tabla ANOVA
Pasos que se realizaron para hallar los valores de la tabla ANOVA:
El cociente F está dado por:
\[ F_0 = \frac{SS_{tratamientos}/(a-1)} {SS_{Error}/(N-a)}= \frac{MS_{tratamientos}} {MS_{Error}}\] donde,
(a-1) son los grados de libertad de los tratamientos, en este caso se realizaron 9 tratamientos, es decir, los grados de libertad son 8.
(N-a) son los grados de libertad del error (72-9=63)
Los grados de libertad totales son N-1, donde N corresponde al número de observaciones en total, en este caso, N= 72 y por tanto, N-1= 71
Para calcular el cuadrado medio \(MS_{tratamientos}\) , se divide la suma de cuadrados entre los grados de libertad. (6000/8 = 750)
- Ahora, para calcular la suma de cuadrados totales, se debe usar la siguiente fórmula:
\[Varianza_{total} = \frac{SC_{totales}} {n-1}\] Se sabe que la varianza total es 883, y los grados de libertad totales son 71, por tanto, la SC totales es 62693.
Sabiendo que la \(SC_{totales}= SC_{tratamientos}+ SC_{Error}\), se puede despejar y hallar la $ SC_{Error}$, lo que sería igual a 56693.
Finalmente, teniendo estos valores se calcula los cuadrados medios del error y el F calculado (0.833)
F_c = 0.833
pf= pf(F_c, 8, 63, lower.tail = F) # Calculo del p_valor
pf # p-valor según el F calculado obtenido ## [1] 0.5770386- Hipótesis nula: Igualdad de los promedios del índice de clorofila en todas las condiciones de tratamiento.
- 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)?
Dado que el \(F_{calculado} < F_{tabulado}\) - \(0.833 < 2.8\) no se rechaza la hipótesis nula de igualdad de promedios, es decir, el promedio de índice de clorofila en los 9 tratamientos es igual. Además, se obtuvo un p_valor de 0.577 (>0.05), confirmando que no se rechaza la hipótesis nula.
- ¿Vale la pena comparar las medias de tratamientos a posteriori del ANOVA (prueba de Tukey)?
No es necesario aplicar la prueba de Tukey para comparar las medias de los tratamientos, puesto que según los datos obtenidos en el ANOVA, no se rechaza la hipótesis nula y por tanto, las medias de todos los tratamientos son iguales.
2. Métodos limpiadores de pelusa en granjas de algodón.
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.
library(readxl)
algodon <- read_excel("D:/Documents/Diseno de experimentos/Datos algodon.xlsx")
algodon## # A tibble: 30 x 3
## Metodo Granjero Pesoperdido
## <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 rowsPara correr posteriormente el ANOVA, se va a convertir el archivo de excel en un data.frame, los métodos y los granjeros se van a convertir en factores y las pérdidas de peso en vector.
algodon = as.data.frame(algodon)
Metodo = factor(algodon$Metodo)
Granjero = factor(algodon$Granjero)
Pesoperdido = as.vector(algodon$Pesoperdido)
str(algodon)## '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 ...
## $ Pesoperdido: num NA 6.75 13.05 10.26 8.01 ...- 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()
- Al momento de correr el ANOVA se usa la funcióm lm, ya que corre los datos como un modelo desbalanceado debido al dato que se perdió del granjero 1 con el método 1.
library(rstatix)## Warning: package 'rstatix' was built under R version 4.0.3##
## Attaching package: 'rstatix'## The following object is masked from 'package:stats':
##
## filtermod_alg =aov(lm(Pesoperdido ~ Granjero+Metodo)) # lm para modelos desbalanceados
alg_anova = Anova(mod_alg, type = 3) # Corre un ANOVA tipo 3
alg_anova## Anova Table (Type III tests)
##
## Response: Pesoperdido
## Sum Sq Df F value Pr(>F)
## (Intercept) 198.62 1 143.8751 2.610e-10 ***
## Granjero 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 ' ' 1Una ventaja de usar el ANOVA tipo III es que éste no depende del tamaño muestral de cada grupo y se puede utilizar para modelos no balanceados, con distinto número de observaciones en cada grupo.
mod_alg2 = lm(Pesoperdido ~ Metodo+Granjero) # Modelo lineal cambiando el orden de los factores
alg_anova2 = Anova(mod_alg2, type = 3) # Se corre un ANOVA tipo 3, cambiando el orden de los factores
alg_anova2## Anova Table (Type III tests)
##
## Response: Pesoperdido
## Sum Sq Df F value Pr(>F)
## (Intercept) 198.62 1 143.8751 2.610e-10 ***
## Metodo 49.12 4 8.8951 0.0003186 ***
## Granjero 139.66 5 20.2331 5.163e-07 ***
## Residuals 26.23 19
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Al cambiar el orden de los efectos en el modelo ANOVA, solamente se intercambia la posición en que aparecen en la tabla del análisis de varianza, sin embargo, se obtienen los mismos valores.
En el ANOVA se observa que existen diferencias entre los métodos de limpieza de las fibras de algodón. Además, no se puede interpretar el p_valor obtenido del granjero, ya que hace referencia al factor de bloqueo.
Sin embargo, se podría cuestionar el uso de las granjas como bloqueo en el estudio, puesto que los seis fardos fueron mezclados a fondo y luego fue procesada una igual cantidad de algodón, usando cada uno de los métodos. Por lo tanto, no se está viendo el impacto de cada una de las granjas, en el proceso.
- 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).
M1_1 = (6.75 + 13.05 + 10.26 + 8.01 + 8.42)/5
M1_1## [1] 9.298El dato faltante para el granjero 1 y el método 1 se va a reemplazar por 9.298
Los nuevos datos para el algodón son:
library(readxl)
algodon2 <- read_excel("D:/Documents/Diseno de experimentos/Datos algodon2.xlsx")
algodon2## # A tibble: 30 x 3
## Metodo Granjero Pesoperdido
## <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 rowsalgodon2 = as.data.frame(algodon2)
Metodo2 = factor(algodon2$Metodo)
Granjero2 = factor(algodon2$Granjero)
Pesoperdido2 = as.vector(algodon2$Pesoperdido)
str(algodon2)## '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 ...
## $ Pesoperdido: num 9.3 6.75 13.05 10.26 8.01 ...Ahora, se va a correr el análisis de varianza pero ya no como un modelo desbalanceado.
mod_aov_alg3 = aov(Pesoperdido2 ~ Granjero2 + Metodo2) # ANOVA para el modelo balanceado
anova_alg3 = aov(mod_aov_alg3) # Anova para el modelo, balanceado
summary(anova_alg3)## Df Sum Sq Mean Sq F value Pr(>F)
## Granjero2 5 139.36 27.873 21.015 2.45e-07 ***
## Metodo2 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 ' ' 1Al comparar el modelo desbalanceado y el balanceado, se observa que en ambos análisis de varianza se rechaza la hipótesis nula de que el peso pérdido de algodón al utilizar los diferentes métodos de limpiador de fibras de algodón es el mismo, es decir, existen diferencias entre ellos.
plot(Pesoperdido~Metodo)
plot(Pesoperdido2~Metodo2)
medias_algodon <- tapply(Pesoperdido2, list(Metodo2), mean)
medias_algodon## M1 M2 M3 M4 M5
## 9.298000 6.195000 7.021667 6.168333 8.760000desviacion_algodon <- tapply(Pesoperdido2, list(Metodo2), sd)
desviacion_algodon## M1 M2 M3 M4 M5
## 2.187468 2.912338 2.013856 2.646971 2.974787coeficiente_variacion = (100*desviacion_algodon)/medias_algodon
coeficiente_variacion## M1 M2 M3 M4 M5
## 23.52622 47.01110 28.68060 42.91226 33.95876Dado que los coeficientes de variación de todos los métodos es mayor al 20%, significaría que la variación de los datos es peligrosa para tomar decisiones, por lo que se recomendaría reevaluar el diseño y la forma en que se tomaron los datos para poder recomendar un método de limpieza del algodón.
3. Datos de carbono orgánico en el suelo a 5 y 10 cm de profundidad.
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%.
Para este ejercicio: Cédula Ana María Montaño: 1010113408 U= 8 T= 1
set.seed(2022)
CO_5 = runif(50, 3, 3.9) # Genera los datos con distribución uniforme para la capa superior
CO_10 = runif(50, 2, 2.3) # Datos para la capa inferior
xy = expand.grid(x=0:sqrt(100), y=0:sqrt(200)) # Se crea la ventana de observación
Suelo_5 = sort.int(CO_5,partial = 33)
Suelo_10 = sort.int(CO_10, partial = 33)
plot(xy)
Carbono_O = c(Suelo_10, Suelo_5)
plot(xy, col = Carbono_O, pch=16)
grafico = plot(xy,cex= Carbono_O,pch=16)
names(xy)## [1] "x" "y"library(plot3D)## Warning: package 'plot3D' was built under R version 4.0.3library(scatterplot3d)## Warning: package 'scatterplot3d' was built under R version 4.0.3xn = xy$x[1:100]
yn = xy$y[1:100]
scatter3D(xn, yn, Carbono_O, ticktype = "detailed", pch = 18, cex =0.8, clab = c("Carbono Orgánico"), cex = 0.2, zlab = "Carbono_O", bty = "g")
- Compare si se encuentran diferencias en la media de carbono entre capas utilizando un nivel de confianza del 95%.
Se puede aplicar una prueba t-students para muestras pareadas, con un nivel de confianza del 95%.
tablaCO = data.frame(Suelo_5, Suelo_10)
tablaCO## Suelo_5 Suelo_10
## 1 3.090582 2.058505
## 2 3.061178 2.085546
## 3 3.108296 2.064168
## 4 3.153794 2.110736
## 5 3.166257 2.039106
## 6 3.230870 2.119438
## 7 3.066869 2.090674
## 8 3.037778 2.118185
## 9 3.146783 2.022874
## 10 3.123485 2.042442
## 11 3.001676 2.020655
## 12 3.143820 2.069504
## 13 3.130267 2.048449
## 14 3.110260 2.119291
## 15 3.548529 2.017852
## 16 3.467457 2.028792
## 17 3.327418 2.083563
## 18 3.463058 2.062073
## 19 3.535903 2.137590
## 20 3.565159 2.173026
## 21 3.472082 2.123820
## 22 3.513294 2.157219
## 23 3.383902 2.138065
## 24 3.397530 2.161123
## 25 3.333285 2.166953
## 26 3.494890 2.164064
## 27 3.483310 2.168711
## 28 3.489420 2.121668
## 29 3.431057 2.129903
## 30 3.364828 2.144229
## 31 3.421189 2.127873
## 32 3.572212 2.121362
## 33 3.573869 2.178137
## 34 3.696249 2.188570
## 35 3.897889 2.221770
## 36 3.582533 2.189953
## 37 3.761429 2.223109
## 38 3.841573 2.234124
## 39 3.724358 2.193088
## 40 3.688558 2.223738
## 41 3.578019 2.222649
## 42 3.681528 2.242517
## 43 3.584443 2.237116
## 44 3.758553 2.257956
## 45 3.780510 2.260114
## 46 3.734380 2.257756
## 47 3.798435 2.249838
## 48 3.765819 2.297042
## 49 3.695896 2.246042
## 50 3.767387 2.239148prueba_tstudent = t.test(tablaCO$Suelo_5,tablaCO$Suelo_10, alternative = "t", paired = T, conf.level = 0.95)
prueba_tstudent##
## Paired t-test
##
## data: tablaCO$Suelo_5 and tablaCO$Suelo_10
## t = 46.825, df = 49, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 1.250864 1.363045
## sample estimates:
## mean of the differences
## 1.306955ifelse(prueba_tstudent$p.value<0.05, "Rechazo Ho", "No rechazo Ho")## [1] "Rechazo Ho"Se encuentran diferencias en las medidas de carbono orgánico a 5cm y 10cm, por lo que se rechaza la hipótesis nula de igualdad de medias entre ambas medidas de carbono. Se observa que a mayor profundidad, los valores de CO disminuyen.
4. Diseño factorial completo \((3^2)\) en arreglo completamente al azar.
El siguiente diseño se corresponde con un factorial completo (\(3^2\)) 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?
- Respuesta:
Este diseño presenta 2 factores, cada uno con 3 niveles. Los tratamientos son 3^2 = 9.
Cada tratamiento se repite 2 veces.
D <- expand.grid(F1= c(3.25, 3.75, 4.25), F2 = c(4,5,6)) # Crea el diseño 3^2
D <- rbind(D,D) # Crea la estructura para dos repeticiones por tratamiento
set.seed(2020)
repetic = paste0("r", 1:6)
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
View(D)library(collapsibleTree)
collapsibleTree(D,hierarchy=c("F1","F2")) # Crea el árbol del diseño que muestra las 18 observaciones. # D corresponde a la variable respuesta- Modelo del diseño:
\[y_{ijk}=\mu+\tau_i+\delta_j+(\tau\delta)_{ij}+\epsilon_{ijk\\i:1\cdots 3\\j:1,3\\k:1\cdots 2}\\Condiciones~Laterales~Respectivas\]
donde
\(Y_{ijk}\) Es la \(ijk-esima\) biomasa para el \(i-esimo\) nivel del factor 1 (F1) y el \(j-esimo\) nivel del factor 2 (F2) y la \(k-esima\) repetición.
\(\mu\) Es la media general de biomasa.
\(\tau_i\) Es el efecto del \(i-esimo\) nivel del F1.
\(\delta_j\) Es el efecto del \(j-esimo\) nivel del F2.
\((\tau\delta)_{ij}\) Es la interacción del \(i-esimo\) nivel del F1 con el \(j-esimo\) nivel del F2.
\(\epsilon_{ijk}\) Es el error residual para el \(i-esimo\) nivel del F1, el \(j-esimo\) nivel del F2 y la \(k-esima\) repetición.
F1: Dosis de un insecticida que se sospecha tiene un efecto de disminución del crecimiento (biomasa)
F2: Número de aplicaciones durante el desarrollo del cultivo.
- Hipótesis:
Primera hipótesis: El efecto de la interacción entre los dos factores, es nulo para todo i,j.
\[H_{01}: (\tau\delta)_{ij}=0~\\ \forall i,j\] Segunda hipótesis: El efecto del factor 1 (Dosis de insecticida) es nulo.
\[H_{02}: \tau_i=0~\\ i=1\cdots 3\]
Tercera hipótesis: El efecto del factor 2 (Aplicaciones) es nulo.
\[H_{03}: \delta_j=0~\\ j= 1,3\]
Ahora, un pequeño análisis descriptivo:
Dosis_insecticida= as.factor(D$F1)
Aplicaciones = as.factor(D$F2)
biomasa_n = as.vector(D$biomasa)
D_NUEVO <- data.frame(Dosis_insecticida, Aplicaciones, biomasa_n)
str(D_NUEVO)## 'data.frame': 18 obs. of 3 variables:
## $ Dosis_insecticida: Factor w/ 3 levels "3.25","3.75",..: 2 1 1 1 1 3 2 2 2 2 ...
## $ Aplicaciones : Factor w/ 3 levels "4","5","6": 1 1 3 2 2 2 3 3 2 2 ...
## $ biomasa_n : num 2.71 2.77 2.14 2.56 2.71 ...lattice::bwplot(D_NUEVO$biomasa_n~D_NUEVO$Aplicaciones|D_NUEVO$Dosis_insecticida) 
Gráficamente, se puede observar que el mayor valor de biomasa podría alcanzarse en la dosis 4.25 con 4 aplicaciones, seguido de la dosis 3.75 con 5 aplicaciones. De forma contraria, el tratamiento que genera menor biomasa es el de dosis 3.25 con 6 aplicaciones.
Si se analiza el patrón de la biomasa dentro de cada dosis, en la de 3.25 la biomasa disminuye a medida que aumentan las aplicaciones. En la de 3.75, la biomasa no sigue un patrón específico y en la de 4.75 la biomasa tiende a disminuir con el aumento de las aplicaciones.
- ANOVA:
aov_mo2 <- aov(biomasa_n~Dosis_insecticida*Aplicaciones) # Corre el análisis de varianza
summary(aov_mo2)## Df Sum Sq Mean Sq F value Pr(>F)
## Dosis_insecticida 2 0.4161 0.20804 1.979 0.194
## Aplicaciones 2 0.3426 0.17128 1.630 0.249
## Dosis_insecticida:Aplicaciones 4 0.1635 0.04087 0.389 0.812
## Residuals 9 0.9459 0.10510Según el análisis de varianza no existe interacción entre las dosis del insecticida y las aplicaciones, puesto que se observa que según el p_valor (0.812>0.05) no se rechaza la hipótesis número 1 que dice que la interacción entre ambos factores es nula.
En cuanto a la segunda y tercera hipótesis, se evidencia que tampoco se rechazan, por lo que el efecto de la dosis del insecticida y de las aplicaciones son nulos, por lo cual, tanto las dosis como las aplicaciones son estadísticamente iguales.
- Prueba de comparación de medias: para confirmar que no hay diferencias entre los tratamientos se realizó la prueba de comparación de medias de Tukey, sin embargo, según el resultado del ANOVA donde las dosis y las aplicaciones son iguales entre sí, no sería necesario correr la prueba Tukey.
prueba_tukey <- aov(biomasa_n~Dosis_insecticida+Aplicaciones)
TukeyHSD(prueba_tukey)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = biomasa_n ~ Dosis_insecticida + Aplicaciones)
##
## $Dosis_insecticida
## 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
##
## $Aplicaciones
## 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.8821270En las comparaciones de medias del efecto de las dosis del insecticida sobre la biomasa se observa que el p-valor ajustado en todos los casos es mayor al 5%, lo que significa que usar una dosis de 3.75 o 3.25 es lo mismo, al igual que en 4.25-3.25 y 4.25-3.75. 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.
Además, en las comparaciones de medias del efecto de las aplicaciones sobre la biomasa, se observa que tienen el mismo comportamiento que las dosis del insecticida y por tanto, su efecto sobre la biomasa es el mismo.
- Ahora, para analizar la interacción calculando las medias, se emplea la función tapply()
t <- tapply(biomasa_n, list(Dosis_insecticida, Aplicaciones), mean) # Calcula las medias
addmargins(t, FUN= mean) # Añades las márgenes de las medias ## 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.940436En el ANOVA se obtuvo que no exite interacción significativa entre los dos factores. Al mirar la tabla de medias se observa que el mejor tratamiento es la dosis 4.25 con 4 aplicaciones, concordando con el mejor tratamiento que muestran las márgenes. Lo anterior, lleva a confirmar una vez más que no existe interacción entre los factores.
Los resultados obtenidos en la tabla de medias, coinciden con lo analizado gráficamente. Por otra parte, al no existir interacción, se puede correr el ANOVA, sin la interacción.
# ANOVA para el modelo sin interacción
anova_sin <- aov(biomasa_n~Dosis_insecticida+Aplicaciones)
summary(anova_sin)## Df Sum Sq Mean Sq F value Pr(>F)
## Dosis_insecticida 2 0.4161 0.20804 2.438 0.126
## Aplicaciones 2 0.3426 0.17128 2.007 0.174
## Residuals 13 1.1094 0.08534Al correr el análisis de varianza sin la interacción, aumentan los grados de libertad de los residuales, pero se llegan a los mismos resultados que el ANOVA anterior.
- Gráfico de interacción con ggplot2.
library(dplyr)##
## 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)
D_NUEVO%>%
group_by(Dosis_insecticida, Aplicaciones) %>%
summarise(medias_groups = mean(biomasa_n)) -> medias2## `summarise()` regrouping output by 'Dosis_insecticida' (override with `.groups` argument)ggplot(medias2)+
aes(x=Dosis_insecticida,y=medias_groups,col = Aplicaciones) +
geom_line(aes(group = Aplicaciones))+
geom_point()
interaction.plot(Dosis_insecticida, Aplicaciones, biomasa_n, main = "Interacción entre Dosis y Aplicaciones", xlab = "Dosis de insecticida", ylab = "biomasa", col = c(1:3))
Según el gráfico obtenido con la librería ggplot2, si existe interacción entre ambos factores, contradiciendo lo obtenido en el ANOVA. Además, se confirma que el mejor tratamiento es el de la dosis 4.25 con 4 aplicaciones.
Sin embargo, es más confiable el resultado obtenido por el ANOVA (posteriormente se revisarán todos los supuestos, validándose el uso de los p_valores del análisis de varianza), ya que en el gráfico se observa interacción, la cual puede deberse a que las medias de todos los tratamientos son muy similares.
CONCLUSIÓN: Si se tienen en cuenta los resultados obtenidos en el ANOVA, agronómicamente se recomiendaría usar la dosis 3.25 con 4 aplicaciones, ya que sería la más eficiente económicamente debido a que en el efecto de los otros tratamientos no se encontraron diferencias significativas, al realizar el ANOVA y la prueba de Tukey para comparar los efectos. Sin emmbargo, de los 9 tratamientos, el que genera mayor biomasa es el de dosis 4.25 con 4 aplicaciones.
- Covariable: Arcilla expansible
Arcilla <- sort.default(runif(18, 0.20, 0.40), decreasing = TRUE) # Crea los datos de la covariable arcilla, con distribución uniforme.
Arc = data.frame(D, Arcilla)
Arc## F1 F2 biomasa Arcilla
## 2 3.75 4 2.708826 0.3985124
## 10 3.25 4 2.772692 0.3768873
## 16 3.25 6 2.143359 0.3726635
## 4 3.25 5 2.560519 0.3345261
## 13 3.25 5 2.708666 0.3302072
## 6 4.25 5 2.773705 0.3228664
## 17 3.75 6 2.770350 0.3227439
## 8 3.75 6 2.832470 0.3153687
## 14 3.75 5 2.898280 0.3111702
## 5 3.75 5 3.359619 0.3060593
## 9 4.25 6 3.054099 0.2940596
## 1 3.25 4 3.157896 0.2848798
## 12 4.25 4 3.487669 0.2775186
## 11 3.75 4 3.451547 0.2590587
## 15 4.25 5 3.016111 0.2575276
## 7 3.25 6 3.042156 0.2457057
## 3 4.25 4 3.200552 0.2334734
## 18 4.25 6 2.989329 0.2114003- Nuevo modelo del diseño
\[y_{ijk}=\mu+\tau_i+\delta_j+(\tau\delta)_{ij}+\beta(x_{ijk}-\bar{x}) +\epsilon_{ijk\\i:1\cdots 3\\j:1,3\\k:1\cdots 2}\\Condiciones~Laterales~Respectivas\]
donde
\(Y_{ijk}\) Es la \(ijk-esima\) observación de biomasa para el \(i-esimo\) nivel del factor 1 (F1) y el \(j-esimo\) nivel del factor 2 (F2) y la \(k-esima\) repetición.
\(\mu\) Es la media general de biomasa.
\(\tau_i\) Es el efecto del \(i-esimo\) nivel del F1.
\(\delta_j\) Es el efecto del \(j-esimo\) nivel del F2.
\((\tau\delta)_{ij}\) Es la interacción del \(i-esimo\) nivel del F1 con el \(j-esimo\) nivel del F2.
\(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}\)
\(\beta\) 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 F1, el \(j-esimo\) nivel del F2 y la \(k-esima\) repetición.
- ANCOVA:
Dosis_Arc = as.factor(Arc$F1)
Aplicaciones_Arc = as.factor(Arc$F2)
biomasa_Arc = as.vector(Arc$biomasa)ancova_arc <- aov(biomasa_Arc~Arcilla+Dosis_Arc+Aplicaciones_Arc, data = Arc)
summary(ancova_arc)## Df Sum Sq Mean Sq F value Pr(>F)
## Arcilla 1 0.8130 0.8130 22.233 0.000501 ***
## Dosis_Arc 2 0.1884 0.0942 2.577 0.117222
## Aplicaciones_Arc 2 0.4279 0.2139 5.850 0.016846 *
## Residuals 12 0.4388 0.0366
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Según el análisis de covarianza, existe evidencia estadística para afirmar que la Arcilla (covariable) influye en la biomasa. Además, se observa que hay diferencias en la biomasa, causadas por las aplicaciones. Esto se debe a que al correr el ANOVA con una covariable se reducen los grados de libertad del error, aumentando de esta forma la potencia para detectar diferencias causadas por los factores, en este caso, de las aplicaciones.
Según, Minitab las covariables se utilizan generalmente en ANOVA y Diseño de Experimentos. En estos modelos, una covariable es cualquier variable continua, que usualmente no se controla durante la recogida de datos. Incluyendo las covariables, el modelo permite incluir y adaptar las variables que se midieron pero no se aleatorizaron o controlaron en el experimento. Al agregar covariables se puede mejorar considerablemente la exactitud del modelo y se puede afectar significativamente los resultados del análisis final. Incluyendo una covariable en el modelo se puede reducir el error en el modelo para incrementar la potencia de las pruebas de los factores.
Revisión de supuestos del ANCOVA:
- Supuesto de linealidad: Para comprobar el supuesto de linealidad se realiza un test de correlación entre la biomasa y la covariable.
plot(Arc$Arcilla, Arc$biomasa) # Gráfico de relación entre la covariable y la variable biomasa
cor.test(Arc$biomasa, Arc$Arcilla)##
## Pearson's product-moment correlation
##
## data: Arc$biomasa and Arc$Arcilla
## t = -3.5112, df = 16, p-value = 0.002894
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.8612947 -0.2786494
## sample estimates:
## cor
## -0.6596974Se observa que existe una correlación negativa moderable (correlación = -0.6596974) entre la variable respuesta biomasa y la covariable arcilla. Además, como el p_valor es menor al 0.05 (0.002894), existe correlación entre ambas variables.
- Homogeneidad de las pendientes de regresión
Con esta suposición se verifica que no haya una interacción significativa entre la covariable y la variable respuesta. Esto se puede evaluar de la siguiente manera:
Arc %>% anova_test(biomasa ~ Arcilla *(F1+F2)+Arcilla*F1*F2)## Coefficient covariances computed by hccm()## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 Arcilla 1 10 15.173 0.003 * 0.603
## 2 F1 1 10 0.807 0.390 0.075
## 3 F2 1 10 11.658 0.007 * 0.538
## 4 Arcilla:F1 1 10 5.183 0.046 * 0.341
## 5 Arcilla:F2 1 10 0.012 0.914 0.001
## 6 F1:F2 1 10 0.456 0.515 0.044
## 7 Arcilla:F1:F2 1 10 0.402 0.540 0.039Se observa que no existe una interacción entre la covariable y los dos factores a la vez, cumpliéndose el segundo supuesto.
- Normalidad de los residuos
# Fit the model, the covariate goes first
modelo11 <- lm(biomasa_Arc ~ Arcilla + (Dosis_Arc+Aplicaciones_Arc), data = Arc)
# Inspect the model diagnostic metrics
model.metrics <- augment(modelo11)shapiro_test(model.metrics$.resid)## # A tibble: 1 x 3
## variable statistic p.value
## <chr> <dbl> <dbl>
## 1 model.metrics$.resid 0.937 0.258La prueba de Shapiro Wilk no fue significativa (p> 0.05), por lo que podemos asumir la normalidad de los residuos. Se cumple el tercer supuesto.
- Homogeneidad de la varianza
Arc$tratamientos <- interaction(Dosis_Arc, Aplicaciones_Arc)
bartlett.test(ancova_arc$residuals~Arc$tratamientos)##
## Bartlett test of homogeneity of variances
##
## data: ancova_arc$residuals by Arc$tratamientos
## Bartlett's K-squared = 9.2079, df = 8, p-value = 0.3251Se cumple el cuarto supuesto de homogeneidad de la varianza puesto que el p_valor es mayor a 0.05, y por lo tanto, no se rechaza la hipótesis nula de que las varianzas son iguales.
- Valores atípicos
model.metrics %>%
filter(abs(.std.resid) > 3) %>%
as.data.frame()## [1] .rownames biomasa_Arc Arcilla Dosis_Arc
## [5] Aplicaciones_Arc .fitted .resid .std.resid
## [9] .hat .sigma .cooksd
## <0 rows> (or 0-length row.names)Se cumple el último supuesto debido a que no existen valores atípicos.
Revisión de supuestos del ANOVA:
- Prueba de normalidad del vector de residuales del modelo
Extraer el vector de residuales
resid= aov_mo2$residuals
mean(resid)## [1] -1.424785e-17var(resid)## [1] 0.0556437sum(resid^2)## [1] 0.9459429shapiro.test(resid)##
## Shapiro-Wilk normality test
##
## data: resid
## W = 0.99029, p-value = 0.999# Claramente se observa que los residuales son normales (p-value > 0.05)Se cumple el primer supuesto, ya que el p_ valor es mayor al 5% (0.999>0.05) y por tanto, los residuales tienen distribución normal.
plot(aov_mo2)
- Prueba de igualdad de varianzas. Homogeneidad de varianzas
tratamientos <- interaction(Dosis_insecticida, Aplicaciones)
bartlett.test(aov_mo2$residuals~tratamientos)##
## Bartlett test of homogeneity of variances
##
## data: aov_mo2$residuals by tratamientos
## Bartlett's K-squared = 7.5942, df = 8, p-value = 0.4741# Claramente se observa que los residuales tienen la misma variación . Permiso para poder usar el p-valor.
# Residuales tienen masomenos la misma variaciónSe cumple el segundo supuesto puesto que el p_valor es mayor a 0.05, y por tanto, no se rechaza la hipótesis nula de que las varianzas son iguales.
- Independencia de los residuales
plot(aov_mo2$residuals, pch=16) # Señal de que los residuales son independientes, no tienen dependencia espacial ni temporal. No siguen ningún patrón. 
Se cumple el tercer supuesto, el de independencia de los residuales.
Cuando las tres condiciones se cumplen, existe una confianza para usar el p_valor obtenido en el ANOVA, por tanto, se verifica que no existe interacción entre las dosis de insecticida y el número de aplicaciones.
CONCLUSIÓN: Al analizar la correlación entre Arcilla (covariable) y la biomasa se observa que si se relacionan de forma lineal y por tanto, si se justifica el uso de la covariable. El uso de la covariable se justifica ya que ayuda a detectar mejor las diferencias entre los tratamientos, en este caso, causadas por las aplicaciones. Al correr el ANOVA no se identificó efecto ni de las dosis ni de las aplicaciones, mientras que cuando se corre el análisis de covarianza, se detecta que existen diferencias entre las aplicaciones.
- Nuevo gráfico de interacción teniendo en cuenta la covariable Arcilla:
Dosis_insecticida2= as.factor(Arc$F1)
Aplicaciones2 = as.factor(Arc$F2)
biomasa_n2 = as.vector(Arc$biomasa)
D_NUEVO2 <- data.frame(Dosis_insecticida2, Aplicaciones2, biomasa_n2)
str(D_NUEVO2)## '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 ...
## $ biomasa_n2 : num 2.71 2.77 2.14 2.56 2.71 ...library(dplyr)
library(ggplot2)
D_NUEVO2 %>%
group_by(Dosis_insecticida2, Aplicaciones2) %>%
summarise(medias_groups = mean(biomasa_n2)) -> medias_nuevas## `summarise()` regrouping output by 'Dosis_insecticida2' (override with `.groups` argument)ggplot(medias_nuevas)+
aes(x=Dosis_insecticida2,y=medias_groups,col = Aplicaciones2) +
geom_line(aes(group = Aplicaciones2))+
geom_point()
No se observan diferencias entre los gráficos de interacción.
5. Diseño anidado escalonado (staggered nested design).
Corresponde a un diseño anidado, factorial incompleto.
a: niveles del factor superior. Finca A y B. a: 20 y a-1: 19
En este ejemplo, el diseño anidado escalonado se utiliza para estimar las fuentes de variabilidad estimando la varianza atribuible a las parcelas, a las parcelas anidadas en las fincas, y a la variedad dentro de la finca.
Se reconoció que este diseño anidado escalonado corresponde a uno de 4 etapas, donde a cada uno de ellas se le asigna una letra de la A a la D.
La respuesta puede ser el rendimiento de la variedad de papa.
library(daewr)## Warning: package 'daewr' was built under R version 4.0.3## Registered S3 method overwritten by 'DoE.base':
## method from
## factorize.factor conf.designdata("polymer")
polymer## lot box prep test strength
## 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
## 81 21 1 1 1 7.60
## 82 21 1 1 2 6.72
## 83 21 1 2 1 6.51
## 84 21 2 1 1 6.35
## 85 22 1 1 1 5.12
## 86 22 1 1 2 5.85
## 87 22 1 2 1 6.31
## 88 22 2 1 1 8.74
## 89 23 1 1 1 5.28
## 90 23 1 1 2 5.73
## 91 23 1 2 1 4.53
## 92 23 2 1 1 5.07
## 93 24 1 1 1 5.44
## 94 24 1 1 2 5.38
## 95 24 1 2 1 4.35
## 96 24 2 1 1 7.04
## 97 25 1 1 1 3.50
## 98 25 1 1 2 3.88
## 99 25 1 2 1 2.57
## 100 25 2 1 1 3.76
## 101 26 1 1 1 4.80
## 102 26 1 1 2 4.46
## 103 26 1 2 1 3.48
## 104 26 2 1 1 3.18
## 105 27 1 1 1 5.35
## 106 27 1 1 2 6.39
## 107 27 1 2 1 4.38
## 108 27 2 1 1 5.50
## 109 28 1 1 1 3.09
## 110 28 1 1 2 3.19
## 111 28 1 2 1 3.79
## 112 28 2 1 1 2.59
## 113 29 1 1 1 5.30
## 114 29 1 1 2 4.72
## 115 29 1 2 1 4.39
## 116 29 2 1 1 6.13
## 117 30 1 1 1 7.09
## 118 30 1 1 2 7.82
## 119 30 1 2 1 5.96
## 120 30 2 1 1 7.14Se creó el data.frame con los primeros 80 datos, cambiando el nombre de las variables
Parcelas =gl(20,4,80)
Finca = rep(c(1, 1, 1, 2), times = 20)
Variedad = c(1,1,2,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2,1)
test = c(1,2,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2,1,1)
Rendimiento = c(9.76, 9.24, 11.91, 9.02, 10.65, 7.77, 10.00, 13.69, 6.50, 6.26, 8.02, 7.95, 8.08, 5.28, 9.15,7.46, 7.84, 5.91, 7.43, 6.11, 9.00, 8.38, 7.01, 8.58, 12.81, 13.58, 11.13, 10.00, 10.62, 11.71, 14.07, 14.56, 4.88, 4.96, 4.08, 4.76, 9.38, 8.02, 6.73, 6.99, 5.91, 5.79, 6.59, 6.55, 7.19, 7.22, 5.77, 8.33, 7.93, 6.48, 8.12, 7.43, 3.70, 2.86, 3.95, 5.92, 4.64, 5.60, 5.96, 5.88, 5.94, 6.28, 4.18, 5.24, 9.50, 8.00, 11.25, 11.14, 10.93, 12.16, 9.51, 12.71, 11.95, 10.58, 16.79, 13.08, 4.34, 5.45, 7.51, 5.21)
P <- data.frame(Parcelas, Finca, Variedad, test, Rendimiento)
P## Parcelas 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.60
## 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# Crea el árbol para el diseño anidado escalonado.
library(collapsibleTree)
collapsibleTree(P,hierarchy=c("Finca","Variedad","test"))Se van a estimar los componentes de la varianza usando el método REML, para ello se usa la función lmer del paquete lm4. Los resultados se muestran a continuación:
library(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 car##
## Attaching package: 'lme4'## The following object is masked from 'package:daewr':
##
## cakemodelo3 <- lmer(Rendimiento ~ 1 + (1|Parcelas) + (1|Parcelas:Finca)+
(1|Parcelas:Finca:Variedad), data = P)## boundary (singular) fit: see ?isSingularsummary(modelo3)## Linear mixed model fit by REML ['lmerMod']
## Formula: Rendimiento ~ 1 + (1 | Parcelas) + (1 | Parcelas:Finca) + (1 |
## Parcelas:Finca:Variedad)
## Data: P
##
## REML criterion at convergence: 325.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.93045 -0.39821 0.01024 0.44045 1.64817
##
## Random effects:
## Groups Name Variance Std.Dev.
## Parcelas:Finca:Variedad (Intercept) 1.2361 1.112
## Parcelas:Finca (Intercept) 0.0000 0.000
## Parcelas (Intercept) 7.0186 2.649
## Residual 0.8743 0.935
## Number of obs: 80, groups:
## Parcelas:Finca:Variedad, 60; Parcelas:Finca, 40; Parcelas, 20
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 8.2362 0.6191 13.3
## convergence code: 0
## boundary (singular) fit: see ?isSingularEl estimador REML para el componente de varianza para la finca dentro de la parcela es igual a 0, por lo tanto, no hay varianza atribuible a las parcelas anidadas en las fincas.
#Porcentaje de la variación de las Parcelas
Vp= (100*7.0186)/(7.0186+1.2361+0.8743)
Vp## [1] 76.88246#Porcentaje de la variación de las variedades dentro de la finca
Vv= (100*1.2361)/(7.0186+1.2361+0.8743)
Vv## [1] 13.54037En los resultados se observa que el 76,88% de la variación total se debe a la variabilidad entre parcelas, mientras que el 13,54% se debe a la variabilidad de las variedades dentro de la finca y la variabilidad dentro de parcela o de finca a finca es insignificante, en este caso 0. Por lo tanto, si los agricultores desean disminuir la variabilidad en el rendimiento de la papa, se deben enfocar en reducir la variación de los factores que cambian entre parcelas, por ejemplo: propiedades físicas y químicas del suelo.
6. Datos de potasio en 8 laboratorios diferentes.
library(asbio)## Warning: package 'asbio' was built under R version 4.0.3## Loading required package: tcltkdata(K)
K## 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 JAnálisis descriptivo
lattice::bwplot(K$K~K$lab)
En este gráfico se observa que el laboratorio que registra los mayores valores de potasio(k) en el suelo fue el E, sin embargo, este a su vez es el que presenta, aparentemente, mayor variabilidad entre los datos. Mientras que los valores del laboratorio G son los más homogéneos. Además, es sospechoso que los resultados obtenidos por el laboratorio H, presenten valores de K por debajo de los registrados por los demás laboratorios.
Además, se observa que los laboratorios B, D y H muestran un dato atípico.
Las muestras, dado que se tomaron de una mezcla de suelo de la finca, deberían presentar valores de potasio similares pero con ciertas diferencias. Sin embargo, los valores obtenidos por cada laboratorio pueden variar debido a factores como: tipo de reactivos usados en el laboratorio, experiencia del laboratorista y calidad de los equipos, asumiendo que la mezcla se realizó de forma homogénea para entregarse a cada laboratorio.
Ahora, se va a realizar un gráfico de densidad ya que es útil para observar la forma de la distribución de la variable potasio.
#Gráfico de densidad
library(ggplot2)
ggplot(data=K, aes(x= K, fill=lab)) +
geom_density(alpha=0.3)
En la gráfica de densidad se visualiza la distribución de los datos en el intervalo de valores de potasio. Se observa que el laboratorio E presenta menor densidad, porque sus datos no están concentrados en un intervalo pequeño de valores de potasio, sino que se distribuyen desde valores iguales a 200 hasta valores de 400. Mientras que el laboratorio G presenta los datos más concentrados, en un rango de valores de potasio (275-320).
medias <- tapply(K$K, list(K$lab), FUN = mean)
medias## B D E F G H I J
## 326.1111 321.1111 316.5556 315.6667 304.1111 229.3333 313.2222 336.2222Para complementar el cálculo de la media, se puede realizar el cálculo del coeficiente de variación para verificar si en verdad el laboratorio J presenta los mayores valores de potasio.
des = tapply(K$K, list(K$lab), sd)
des## B D E F G H I J
## 44.40564 29.25510 80.35097 25.04995 15.37404 23.89037 31.39577 21.60890Se observa que aunque el laboratorio J no es el que presenta menor desviación, si es uno de los que presenta valores más homogéneos, lo que indica más confiabilidad a los resultados de este laboratorio. Además, se confirma lo que se interpretó en la gráfica de cajas, donde el laboratorio E presenta la mayor desviación, y como se observa en el gráfico de densidad, es el más disperso.
coef_var = 100 * des/medias
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.426969El coeficiente de variación nos ayuda a descartar los posibles mejores promedios como el caso del E, pero que quizá tenga mucha dispersión de los datos. Si tomamos la regla de que un CV mayor al 20% puede resultar peligroso para la toma de decisiones, observamos que solo el laboratorio E supera este valor y se confirma que el laboratorio J es el mejor, ya que presentó el segundo coeficiente de variación más bajo y el mayor promedio, seguido del B y el D.
- Enfoque no paramétrico: Prueba de Krustall-Wallis - Análisis de varianza de rangos.
En esta prueba se van a comparar los promedios de los rangos para cada laboratorio.
\[H_0: \mu_{rangoB}=\mu_{rangoD}=\mu_{rangoE}=\mu_{rangoF}=\mu_{rangoG}=\mu_{rangoH}=\mu_{rangoI}=\mu_{rangoJ}\]
K$K## [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 353which.min(K$K)## [1] 24K$rangos <-rank(K$K)
tapply(K$rangos, K$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.944444Se observa que los rangos promedio son muy similares en los laboratorio B,D y E.
boxplot(K$rangos~K$lab)
Gráficamente, se observa que hay diferencias entre los rangos de cada laboratorio, sin embargo los valores de los B,D,E,F son muy similares, mientras que el H si difiere en mayor grado.
kruskal.test(K$K~K$lab) # Corre la prueba de Kruskal-Wallis##
## Kruskal-Wallis rank sum test
##
## data: K$K by K$lab
## Kruskal-Wallis chi-squared = 24.482, df = 7, p-value = 0.000937Con el p_valor obtenido de 0.000937 (<0.05) se rechaza la hipótesis nula de igualdad de rangos entre los laboratorios, y se confirma lo que se interpretó gráficamente, de que los rangos de los promedios asignados a cada uno de los 8 laboratorios son diferentes.
Ahora, con la prueba de Nemenyi se pueden comparar las medias por pares de laboratorio.
PMCMR::posthoc.kruskal.nemenyi.test(K$K~K$lab) # Corre la 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: K$K by K$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 la prueba de Nemenyi se observa mejor si existen diferencias entre cada par de laboratorios. Se observa que en la mayoría de parejas se obtiene un p-valor mayor al 5%. Sin embargo, en otros casos, los laboratorios si presentan diferencias, principalmente con el laboratorio H que al compararse con todos, presentó un p-valor menor de 0.05. Por lo tanto, se rechazasss la hipótesis nula que dice que todos los laboratorios son iguales.
7. Strip plot design. Diseño de parcelas divididas en bloques completos.
En este ejemplo, se desea evaluar el peso del cultivo de tomate de un ensayo de campo en parcelas divididas utilizando 3 variedades de tomate y 4 niveles de fertilización con Fósforo. El experimento se dispuso en 4 bloques de 3 parcelas principales, cada una de ellas dividida en 4 subparcelas. Las variedades fueron asignadas a las parcelas principales y los tratamientos de fertilización con fósforo a las subparcelas.
En este caso, los investigadores tienen interés de saber cuál tratamiento de fósforo le otorga más peso a los frutos del tomate, por eso, se colocan las dosis de fósforo como subparcelas, ya que los grados de libertad de los residuales del fósforo serán mayores que en la parcela principal (variedad) y por tanto, va a ser más potente detectar diferencias entre estos tratamientos. En consecuencia, se colocan las variedades en la parcela principal.
Además, las variedades que se estudiaron presentan características similares en cuanto al peso del fruto.
set.seed(2025)
Bloqueo = gl(4, 12, 48)
variedad = gl(3, 4,48, c("V1","V2","V3"))
Fosforo = gl(4,1,48, c(0.0, 0.3, 0.6, 0.9))
Peso = rnorm(48, 110, 0.5)
tomate = data.frame(Bloqueo, variedad, Fosforo, Peso)
tomate## Bloqueo variedad Fosforo Peso
## 1 1 V1 0 110.3104
## 2 1 V1 0.3 110.0178
## 3 1 V1 0.6 110.3866
## 4 1 V1 0.9 110.6362
## 5 1 V2 0 110.1855
## 6 1 V2 0.3 109.9186
## 7 1 V2 0.6 110.1986
## 8 1 V2 0.9 109.9600
## 9 1 V3 0 109.8275
## 10 1 V3 0.3 110.3511
## 11 1 V3 0.6 109.8022
## 12 1 V3 0.9 109.1225
## 13 2 V1 0 109.7895
## 14 2 V1 0.3 110.3825
## 15 2 V1 0.6 110.5331
## 16 2 V1 0.9 109.9032
## 17 2 V2 0 109.9739
## 18 2 V2 0.3 110.2088
## 19 2 V2 0.6 109.3437
## 20 2 V2 0.9 111.2145
## 21 2 V3 0 109.8777
## 22 2 V3 0.3 110.3669
## 23 2 V3 0.6 111.4264
## 24 2 V3 0.9 110.3358
## 25 3 V1 0 109.9421
## 26 3 V1 0.3 109.3493
## 27 3 V1 0.6 109.5411
## 28 3 V1 0.9 110.0254
## 29 3 V2 0 109.8099
## 30 3 V2 0.3 110.3613
## 31 3 V2 0.6 109.4871
## 32 3 V2 0.9 110.9099
## 33 3 V3 0 109.3648
## 34 3 V3 0.3 110.2944
## 35 3 V3 0.6 109.6829
## 36 3 V3 0.9 109.3671
## 37 4 V1 0 110.1419
## 38 4 V1 0.3 109.2591
## 39 4 V1 0.6 110.4559
## 40 4 V1 0.9 110.4857
## 41 4 V2 0 110.6366
## 42 4 V2 0.3 110.7311
## 43 4 V2 0.6 109.6945
## 44 4 V2 0.9 109.8197
## 45 4 V3 0 109.2354
## 46 4 V3 0.3 110.5666
## 47 4 V3 0.6 109.7796
## 48 4 V3 0.9 110.0648- Modelo del diseño (strip-plot)
\[y_{ijk}=\mu+\alpha_i+\gamma_k+\eta_{ik}+\beta_j+\ (\alpha\beta)_{ij} +\epsilon_{ijk\\i:1\cdots 3\\j:1,4\\k:1\cdots 4}\\Condiciones~Laterales~Respectivas\]
donde
\(Y_{ijk}\) es el peso del fruto de tomate
\(\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 del fósforo
\((\alpha\beta)_{ij}\) es el efecto de la interacción entre ambos factores
\(\epsilon_{ijk}\) es el error en la subparcela
str(tomate)## 'data.frame': 48 obs. of 4 variables:
## $ Bloqueo : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
## $ variedad: Factor w/ 3 levels "V1","V2","V3": 1 1 1 1 2 2 2 2 3 3 ...
## $ Fosforo : Factor w/ 4 levels "0","0.3","0.6",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ Peso : num 110 110 110 111 110 ...variedad2 = as.factor(tomate$variedad)
Fosforo2 = as.factor(tomate$Fosforo)
Bloqueo2 = as.numeric(tomate$Bloqueo)
str(tomate)## 'data.frame': 48 obs. of 4 variables:
## $ Bloqueo : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
## $ variedad: Factor w/ 3 levels "V1","V2","V3": 1 1 1 1 2 2 2 2 3 3 ...
## $ Fosforo : Factor w/ 4 levels "0","0.3","0.6",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ Peso : num 110 110 110 111 110 ...ANOVA- con la librería agricolae
attach(tomate)
library(agricolae)
model = with(tomate, strip.plot(Bloqueo2, variedad2, Fosforo2, Peso))##
## ANALYSIS STRIP PLOT: Peso
## Class level information
##
## variedad2 : V1 V2 V3
## Fosforo2 : 0 0.3 0.6 0.9
## Bloqueo2 : 1 2 3 4
##
## Number of observations: 48
##
## model Y: Peso ~ Bloqueo2 + variedad2 + Ea + Fosforo2 + Eb + Fosforo2:variedad2 + Ec
##
## Analysis of Variance Table
##
## Response: Peso
## Df Sum Sq Mean Sq F value Pr(>F)
## Bloqueo2 3 1.1368 0.37893 1.4925 0.25033
## variedad2 2 0.2807 0.14035 0.6126 0.57267
## Ea 6 1.3746 0.22909 0.9024 0.51463
## Fosforo2 3 0.4357 0.14525 1.3300 0.32436
## Eb 9 0.9829 0.10921 0.4302 0.90140
## Fosforo2:variedad2 6 3.3721 0.56202 2.2137 0.08943 .
## Ec 18 4.5699 0.25388
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## cv(a) = 0.4 %, cv(b) = 0.3 %, cv(c) = 0.5 %, Mean = 110.0641Según el ANOVA obtenido con la función strip.plot de la librería agricolae se observa que no existe interacción entre ambos factores (dosis de fósforo y las variedades). Asimismo, ni las dosis de fósforo ni la variedad tienen un efecto en el peso del tomate.
Según el ANOVA, no se justificaría el uso de bloqueo ya que no presentó un efecto significativo en el peso del tomate. Sin embargo, se debe considerar inicialmente para evitar tomar decisiones erróneas.
- Añadir covariable: En este caso, se escogió como covariable la salinidad del suelo, medida como conductividad eléctrica (dS/m). Debido a que la solución salina reduce marcadamente el peso seco de la planta, así como el tamaño del fruto y altura.
Salinidad <- sort.default(runif(48, 3, 6), decreasing = TRUE) # Crea los datos de la covariable salinidad , con distribución uniforme.
S <- data.frame(tomate, Salinidad)
S## Bloqueo variedad Fosforo Peso Salinidad
## 1 1 V1 0 110.3104 5.889013
## 2 1 V1 0.3 110.0178 5.880791
## 3 1 V1 0.6 110.3866 5.870994
## 4 1 V1 0.9 110.6362 5.789767
## 5 1 V2 0 110.1855 5.753410
## 6 1 V2 0.3 109.9186 5.745401
## 7 1 V2 0.6 110.1986 5.697654
## 8 1 V2 0.9 109.9600 5.585155
## 9 1 V3 0 109.8275 5.571172
## 10 1 V3 0.3 110.3511 5.426116
## 11 1 V3 0.6 109.8022 5.424559
## 12 1 V3 0.9 109.1225 5.422066
## 13 2 V1 0 109.7895 5.412673
## 14 2 V1 0.3 110.3825 5.376925
## 15 2 V1 0.6 110.5331 5.215864
## 16 2 V1 0.9 109.9032 5.205366
## 17 2 V2 0 109.9739 5.171802
## 18 2 V2 0.3 110.2088 5.166686
## 19 2 V2 0.6 109.3437 5.161563
## 20 2 V2 0.9 111.2145 5.151789
## 21 2 V3 0 109.8777 5.147222
## 22 2 V3 0.3 110.3669 5.080204
## 23 2 V3 0.6 111.4264 5.052542
## 24 2 V3 0.9 110.3358 5.007446
## 25 3 V1 0 109.9421 4.996553
## 26 3 V1 0.3 109.3493 4.933923
## 27 3 V1 0.6 109.5411 4.920833
## 28 3 V1 0.9 110.0254 4.824139
## 29 3 V2 0 109.8099 4.660071
## 30 3 V2 0.3 110.3613 4.446412
## 31 3 V2 0.6 109.4871 4.390305
## 32 3 V2 0.9 110.9099 4.306982
## 33 3 V3 0 109.3648 4.306725
## 34 3 V3 0.3 110.2944 4.301333
## 35 3 V3 0.6 109.6829 4.036574
## 36 3 V3 0.9 109.3671 4.021880
## 37 4 V1 0 110.1419 3.810599
## 38 4 V1 0.3 109.2591 3.733043
## 39 4 V1 0.6 110.4559 3.683046
## 40 4 V1 0.9 110.4857 3.608470
## 41 4 V2 0 110.6366 3.550137
## 42 4 V2 0.3 110.7311 3.498575
## 43 4 V2 0.6 109.6945 3.408977
## 44 4 V2 0.9 109.8197 3.402388
## 45 4 V3 0 109.2354 3.376648
## 46 4 V3 0.3 110.5666 3.360496
## 47 4 V3 0.6 109.7796 3.329165
## 48 4 V3 0.9 110.0648 3.021387- ANCOVA: Análisis de covarianza
plot(S$Salinidad, S$Peso)
cor.test(S$Salinidad, S$Peso)##
## Pearson's product-moment correlation
##
## data: S$Salinidad and S$Peso
## t = 0.57778, df = 46, p-value = 0.5662
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.2041775 0.3603259
## sample estimates:
## cor
## 0.08488151En el test de correlación se observa que la covariable salinidad y la variable respuesta peso del tomate, no tienen una correlación significativa. Además, al obtener un p_valor (0.5662) mayor al 5% se confirma que no existe correlación entre ambas variables.
Al no existir dicha correlación no es necesario correr el análisis de covarianza, ya que los resultados no serían significativos porque la salinidad no influiria en el peso del tomate.
Además, no se justifica utilizar la covariable, ya que no existe correlación entre la salinidad y el peso del tomate.
- Nuevo modelo sin interacción:
\[y_{ijk}=\mu+\alpha_i+\gamma_k+\eta_{ik}+\beta_j+\epsilon_{ijk\\i:1\cdots 3\\j:1,4\\k:1\cdots 4}\\Condiciones~Laterales~Respectivas\]
donde
\(Y_{ijk}\) es el peso del fruto de tomate
\(\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 del fósforo
\(\epsilon_{ijk}\) es el error en la subparcela
Se escribe nuevamente el modelo del diseño, pero quitando la interacción entre ambos factores, ya que según los resultados del análisis de varianza, esta interacción no existe. Además, no es posible escribir el modelo con los términos que presenten un p_valor menor al 6%, ya que al no existir efectos significativos en el ANOVA, todos los p_valores son mayores a 0.06.
- Gráfico:
library(lattice)
with(tomate, xyplot(Peso ~ Fosforo2|variedad2, groups = Bloqueo2, aspect = "xy", type = "o"))
Aspectos a tener en cuenta: Se distingue cada bloque teniendo en cuenta los datos creados en el data.frame.
El bloque 1: corresponde al color azul
El bloque 2: corresponde al color fucsia
El bloque 3: corresponde al color verde
El bloque 4: corresponde al color rojo
Gráficamente, se observa que el mejor tratamiento corresponde al bloque 2 (línea fucsia), para la variedad 3 y una dosis de fósforo de 0.6, ya que en este tratamiento se alcanza el mayor peso para el tomate, presentando las mejores condiciones agronómicamente.
Dado que la covariable no presenta un papel significativo, no se realizó el análisis gráfico respectivo.
8. Resumen
Artículo 1: El uso de los diseños en parcelas divididas.
Un experimento de parcela dividida es una buena opción cuando se tiene una situación donde alguno de los factores tiene niveles que son difíciles de aleatorizar o cambiar. Los factores que son difíciles de aleatorizar se le designan a la parcela principal, mientras que los más fáciles de aleatorizar se designan como factores de subparcela.
Este diseño se originó en un entorno agrícola, de allí que los términos empleados en su análisis hagan referencia a parcelas de tierra. Sin embargo, también se emplean en otros ámbitos como el industrial, donde, por ejemplo, los factores difíciles de cambiar pueden ser la temperatura de un horno o la configuración de un equipo especializado.
En cuanto al número de parcelas completas este corresponderá al número de veces que se repita la aplicación del factor de parcela completa, mientras que el número de subparcelas va a ser equivalente al número total de experimentos.
Adicionalmente, al analizar un diseño de parcelas divididas hay que tener cuidado, y no confundirlo con un diseño completamente al azar, ya que las parcelas divididas pueden pasar inadvertidas cuando los datos no estén bien aleatorizados y por tanto, pueden parecer independientes provocando la determinación errónea del efecto de un factor.
Cuando se ejecuta un diseño de parcelas divididas hay que realizar dos aleatorizaciones, una la que determine el orden en el que se ejecutan las parcelas completas y la segunda, dentro de cada parcela completa. Asimismo, en el modelo estadístico del diseño de parcelas divididas, van a existir dos términos de error, uno atribuible a la parcela principal y el otro a la subparcela.
Por otra parte, aunque en muchas ocasiones estos diseños se consideran más complejos, pueden llegar a producir estimaciones de parámetros más precisas que si se corre como un diseño completamente al azar.
Igualmente, en término de costos este diseño puede resultar más útil, puesto que, se puede seleccionar el número de parcelas completas que le represente un menor costo al investigador; y este costo se puede calcular como la suma ponderada del número de parcelas completas y el número total de subparcelas.
Además, para elegir un buen diseño se deben analizar tanto características cualitativas (ej.el número de niveles de cada factor), como características cuantitativas (ej. términos de error de la subparcela). Finalmente, el investigador debe definir cuál es el óptimo para su experimento así como los criterios más importantes para alcanzar los objetivos de su proyecto.
Artículo 2: Unidad experimental y unidad de observación.
Al momento de plantear un diseño experimental se debe definir de forma adecuada la unidad experimental, la cual no siempre resulta ser tan clara. Para ello, se debe tener en cuenta que la unidad experimental es aquella que recibe el tratamiento. También se conoce como unidad de replicación y corresponde a la entidad más pequeñas que se destina a un tratamiento en particular y que es independiente a todas las demás unidades.
También, se define una réplica del experimento como cada vez que se aplica todo el conjunto de niveles del tratamiento.
Por otro lado, muchas veces se puede confundir la unidad experimental con la unidad de muestreo, lo que conlleva a obtener datos erróneos en el ANOVA, puesto que al tratar una unidad de muestreo como una experimental, los grados de libertad del error van a aumentar, reduciendo la suma de cuadrados medios y produciendo un estadístico F más grande.
El coeficiente de variación nos ayuda a descartar los posibles mejores promedios como el caso del E, pero que quizá tenga mucha dispersión de los datos. Si tomamos la regla de que un CV mayor al 20% puede resultar peligroso para la toma de decisiones, observamos que solo el laboratorio E supera este valor y se confirma que el laboratorio J es el mejor, ya que presentó el segundo coeficiente de variación más bajo y el mayor promedio, seguido del B y el D.
- Short communication: On recognizing the proper experimental unit in animal studies in the dairy sciences:
Las unidades experimentales generalmente se consideran intercambiables, lo que implica estadísticamente que las unidades no difieren de manera fundamental y por tanto, se obtienen inferencias confiables independientemente del tratamiento que se le asigne a cada unidad.
Puede existir variabilidad a nivel de unidades experimentales y dentro de ellas. Para separar de forma eficaz los efectos del tratamiento de los efectos de la unidad experimental se deberían replicar los tratamientos en más unidades experimentales.
Ahora, se define una unidad de observación o unidad de muestreo como aquella entidad física sobre la cual se mide un resultado de interés de un experimento.
Otro concepto de interés es el de: estructura de datos jerárquica, la cual corresponde a una configuración de los datos donde las observaciones no son independientes entre sí y por tanto, presentan determinada correlación dada por el diseño experimental anidado. Cuando las unidades de observación se anidan dentro de una unidad experimental,las unidades de observación se conocen como submuestras o réplicas técnicas, para indicar que estas observaciones están correlacionadas y, por lo tanto, no constituyen una verdadera replicación independiente.
A continuación, se mencionan 2 razones que respaldan el hecho de realizar una adecuada replicación como requisito para obtener experimentos válidos :
Los resultados son reproducibles, al menos en las condiciones experimentales especificadas.
La variabilidad entre las unidades experimentales que definen el error experimental se estima adecuadamente y, por tanto, las pruebas de hipótesis posteriores se calibran adecuadamente.
Igualmente, existen modelo mixtos que son capaces de reconocer al mismo tiempo múltiples fuentes de variación aleatoria en un conjunto de datos, evaluando así la variabilidad a nivel de unidad de observación y la variabilidad a nivel de unidad experimental; utilizando una u otra como error experimental para un tratamiento en específico. Asimismo, estos modelos reconocen adecuadamente la unidad experimental dentro del diseño de interés, lo que lleva a establecer un adecuado nivel de error experimental y, por tanto, a obtener la prueba de hipótesis adecuada.
Finalmente, dependiendo el interés y el modelo o diseño que formule el investigador, es de gran importancia que reconozca cuál va a ser la unidad experimental y cuál la unidad de observación, o si van a ser las mismas. Además, para establecer de forma adecuada la unidad experimental en un estudio, se debe basar en una pregunta de investigación específica.
3. Artículo: Fundamentos del diseño experimental: pautas para diseñar experimentos exitosos.
Existen cuatro principios básicos del diseño experimental: replicación, aleatorización, bloqueo y tamaño de las unidades experimentales, para resolver problemas reales, pensando siempre que el fracaso no es una opción.
Al diseñar experimentos se relacionan las ciencias biológicas con la matemáticas. Dentro de este campo, se pueden distinguir tres tipos de experimentos:
Experimentos de observación: para medir una constante asumida.
Experimentos de medición: para medir las propiedades de una población.
Experimentos comparativos: se pretende comparar dos o más prácticas, que puedan tener alguna relevancia importante en las investigaciones.
En un experimento comparativo se sigue la teoría de la investigación científica donde es importante seguir estos pasos:
Primero, se realiza la formulación de pregunta o hipótesis. Luego, se plantea un modelo basado en el tema de interés, el cual debe ser llevado a un modelo estadístico, que implica a su vez un diseño estadístico al cual se le va a aplicar el análisis estadístico. Para finalmente, interpretar los resultados y obtener nuevas hipótesis y preguntas. Lo anterior, permite diseñar en un futuro experimentos más eficientes.
Por otra parte, cuando en un diseño se obtiene un p valor de 0.5 puede deberse a:
- Un experimento mal diseñado con poder insuficiente para detectar una diferencia entre dos medias.
2.Tratamientos mal diseñados que no reflejan adecuadamente la pregunta o hipótesis inicial.
Se realizó un experimento sin una supervisión adecuada sobre los protocolos de tratamiento y recopilación de datos.
La falta de verdaderas diferencias entre las medias de los tratamientos.
Pautas a tener en cuenta, para que un expermiento sea más eficiente.
Replicación:
Las funciones de la replicación son: proporcionar un mecanismo para estimar el error experimental, también para aumentar la presición del experimento, aumenta el alcance de la inferencia del experimento y va a afectar el control del error. Además, la replicación debe aplicarse a la unidad experimental y debe ocurrir a una escala espacial y temporal que coincida con la aplicación de los tratamientos.
Asimismo, debe escogerse una escala de replicación adecuada, ya que se debe evitar confundir las diferencias de tratamiento con otros factores que pueden variar entre las unidades experimentales.
- Cuando no se realiza la replicación en una escala adecuada, se conoce como pseudoreplicación. Lo anterior, provoca una confusión creada por la falta de replicación del tratamiento que se va a extender a todo el experimento.
Los principales pasos de la replicación son: definir la unidad experimental, replicar los tratamientos en una escala adecuada, si se desea diseñar niveles adicionales de replicación.
Para saber cuántos recursos se deben usar para lograr la replicación en cada uno de los niveles deseados, se tienen que considerar el número de réplicas, ya que tiene un efecto directo sobre la precisión y la capacidad de detectar diferencias entre tratamientos.
Experimentos no replicados: En estos experimentos se realizan diseños aumentados con gran cantidad de tratamientos que no son replicados. En ocasiones, las replicaciones no se realizan por temas económicos.
Aleatorización: Cuando se elige una muestra de una población, esta debe ser aleatoria y representativa.
La aleatorización proporciona dos aspectos clave en los diseños experimentales: primero, la estimación no sesgada de las medias del tratamiento y los errores experimentales; y una precaución contra una alteración que pueda surgir o no.
Bloqueo: El bloqueo se utiliza en el diseño de experimentos con los siguientes propósitos:
Por precisión, para crear grupos de unidades experimentales más homogéneos.
Por conveniencia, para permitir diferentes tamaños de unidades experimentales.
Un ejemplo de bloqueo es el diseño de bloques completos al azar (RCBD) donde el tamaño del bloque debe ser igual al número de tratamientos.
Tamaño de las unidades experimentales: Las unidades experimentales se pueden considerar relativamente como pequeñas o grandes, pero esto depende de muchos factores como la especie, el tipo de suelo o las mediciones. En las unidades “pequeñas”, cualquier cambio en el tamaño tendrá gran efecto en la varianza media. Contrariamente, si se emplea una unidad “grande” los cambios en el tamaño pueden tener poco o ningún impacto en la varianza de la unidad media.
9. Artículo:
Chlorophyll a fluorescence and development of zucchini plants under nitrogen and silicon fertilization.
Fluorescencia de la clorofila a y desarrollo de plantas de calabacín bajo fertilización con nitrógeno y silicio.
Objetivo: Evaluar la interacción entre la fertilización con nitrógeno y silicio sobre el crecimiento, el índice de clorofila y la fluorescencia de la clorofila a de plantas de calabacín.
- La mención de la estructura factorial
Los tratamientos se distribuyeron en un esquema de parcelas divididas en un diseño de bloques al azar con tres repeticiones y una planta por repetición. La parcela se formó por niveles de silicio (0 y 6 g/planta) y las subparcelas constituidas por cinco niveles de nitrógeno (30, 60, 90, 120 y 150 kg ha), totalizando 30 parcelas experimentales.
Dos factores: Nitrógeno con 5 niveles y Silicio con 2 niveles.
Tratamientos: 10 tratamientos para cada variable respuesta que analizaron (masa seca de la hoja, masa seca del tallo, masa seca total, clorofila a, clorofila b, clorofila total, fluorescencia de clorofila).
- La razón de colocar cada factor en la parcela principal o en la subparcela.
Se analiza que se colocó el silicio en la parcela principal, ya que al aplicarlo mediante microrriego, va a ser más difícil aleatorizar los niveles de silicio. Asimismo, queda como subparcela las dosis de nitrógeno (no se especifica como lo aplican).
- La revisión de supuestos para el análisis de varianza
No hay revisión de los supuestos del ANOVA en el artículo.
- La tabla del análisis de varianza
El análisis de varianza (Tabla 3) está bien, ya que sólo analizaron la interacción entre el Ni y el Si, ya que se obtuvo una interacción estadísticamente significativa entre estos factores, a un nivel de probabilidad de 1%. Sin embargo, no muestran explícitamente los resultados de los cocientes F y p_valores, concluyendo a partir de los cuadrados medios. En cuanto a los grados de libertad, primero en los bloques debería ser 2, ya que se realizaron 3 repeticiones, para el Silicio debe ser 1, ya que son dos niveles de Si.
En el análisis de varianza concluyen a partir de los cuadrados medios obtenidos para la masa seca de la hoja, la masa seca del tallo y la masa seca total, tomando estas variables como respuesta.
La suma de grados de libertad debe ser igual a n-1, en este caso, 30-1= 29, que coincide con la suma de los grados de libertad que aparecen en la tabla de análisis de varianza. En el análisis de varianza del efecto del N y el Si en el contenido de clorofila y los parámetros de fluorescencia, se observa que existe interacción entre ambos factores (N y Si). Sin embargo, en el artículo analizan y afirman que el N tiene efecto significativo sobre el contenido de clorofila, analizando el factor por separada, lo cual es incorrecto, ya que presenta interacción con el Si.
Los datos se sometieron a un análisis de varianza mediante la prueba F con una probabilidad del 5%. Se realizaron análisis de regresión polinomial para el factor nitrógeno y la prueba de Tukey para el factor de silicio en caso de significancia. Se utilizó el programa estadístico R
- El uso de muchos análisis de varianzas en lugar de uno solo multivariante
En el artículo se realizaron dos análisis de varianza, el primero para la masa seca, y el segundo para el contenido de clorofila. Sin embargo, si se realiza un sólo análisis multivariante sería más sencillo concluir que tratamiento es el más adecuado para las dos variables al tiempo.
- El método de comparaciones de medias después del Anova
En cuanto a la comparación de pares de medias, aplicando una prueba de Tukey al 5%, para la masa seca de la hoja se observa que son estadísticamente iguales solamente los tratamientos de 120 N con 0 Si y 120 N y 6 con Si, ya que muestran las mismas letras. Para la masa seca de la hoja, el mejor tratamiento es 150 de N y 0 de Si. Media de 355. Para la masa seca del tallo, ningún par de medias son estadísticamente iguales y se observa que el mejor tratamiento corresponde a 90Kg de N y 0 de Si, ya que se obtiene la mayor masa media (46). Si se observa la masa seca total, esta corresponde al tratamiento de 150 de N con 0 Si (375).
La prueba de Tukey se aplica de manera análoga con el índice de clorofila a, clorofila b, clorofila total y fluorescencia de clorofila.
- La interacción de factores
Según el ANOVA si existe interacción entre el N y el Si.
- La presencia de bloques
Existen 3 bloques.
- El balanceo o desbalanceo.
No lo mencionan, pero se puede considerar como un diseño balanceado porque en el artículo no especifican la presencia de algún dato faltante.
- La definición clara de la unidad experimental.
La unidad experimental se define como aquella a la que se aplican los tratamientos, en este caso, se tienen 30 parcelas experimentales. La unidad experimental sería cada planta que se encuentra en cada una de las 30 parcelas experimentales y que recibe un tratamiento de Si y N.
- Software utilizado y librería específica (en caso de ser R)
Utilizaron R pero no mencionan la librería utilizada.