Parcial - D. de experimentos 2020-2
Noviembre 2020
Integrantes : Paula Natalia Paez Monroy, Angie Alejandra Juyo Buitrago
Nota: Número usado como U :3, número usado como T:0
Taller - Parcial diseño de experimentos
1. En un estudio conducido en ambiente controlado se tuvieron 72 macetas, cada una con una planta a la que a cierta edad se le midió el contenido de clorofila (índice de clorofila) con un sensor (SPAD). El total de macetas se correspondió con 9 tratamientos asociados a estrés hídrico. Se sabe que la varianza de las 72 observaciones es 833. Con esta información complete la tabla del ANOVA que se muestra a continuación:
Tabla 1 Tabla de analisis de varianza para el primer punto.
A.Si el F tabulado es 2.8. ¿qué puede decirse acerca de la Hipótesis nula de igualdad de los promedios del índice en todas las condiciones de tratamiento (use el p valor sí como el cociente F calculado de la tabla para concluir)?
\[H_o: \mu_1 = \mu_2 = \mu_3 ... = \mu_9\] \[H_a: H_o \ es \ falsa\] RTA: Al tener un F tabulado mayor que el calculado, permite afirmar que los nueve tratamientos se comportan estadísticamente iguales, ya que el F calculado cae en la zona de no rechazo, por lo cual se aceptaría la hipotesis nula, que dicta que los trantamientos asociados a estres hidrico no son estadisticamente diferentes. Sin embargo, el valor de F calculado es 0.889, esto es señal de que lo que esta causando mayor variabilidad son las repeticiones, por lo tanto se puede afirmar que hay problemas en el ensayo realizado. Posiblemente, no existan diferencias en los tratamientos, pero al tener un F calculado de esta magnitud las repeticiones son las causantes de la variacion, y llegar a la conclusion de que los tratamientos son iguales es preocupante.
Calculo del p.valor
pf(q = 0.926, df1 = 8, df2 = 63, lower.tail = F)## [1] 0.5014761Adicionalmente, se puede observar que, el p.valor es mayor al 5%, por lo cual no se rechaza la hipotesis nula, como ya se había mencionado anteriormente, por lo tanto en relación a la pregunta:¿vale la pena comparar las medias de tratamientos a posteriori del ANOVA (prueba de Tukey)? En vista, del F calculado y el p.valor, no hay que realizar prueba de comparacion de medias, ya que los valores se comportan estadísticamente iguales.
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 investigar 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 de algodón. 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.
Figura 1 Datos del experimento para el segundo punto.
library(readxl) # Improtar datos de excel
datos_peso <- read_excel("datos_g.xlsx")df2 = data.frame(datos_peso) # Generar data.frame
head(df2)## Granja Método Pérdida
## 1 1 M2 5.54
## 2 1 M3 7.67
## 3 1 M4 7.89
## 4 1 M5 9.27
## 5 2 M1 6.75
## 6 2 M2 3.53metodo <- as.factor(df2$Método); metodo## [1] M2 M3 M4 M5 M1 M2 M3 M4 M5 M1 M2 M3 M4 M5 M1 M2 M3 M4 M5 M1 M2 M3 M4 M5 M1
## [26] M2 M3 M4 M5
## Levels: M1 M2 M3 M4 M5granja <-as.factor(df2$Granja) ; granja## [1] 1 1 1 1 2 2 2 2 2 3 3 3 3 3 4 4 4 4 4 5 5 5 5 5 6 6 6 6 6
## Levels: 1 2 3 4 5 6A.Realice el ANOVA para este diseño recordando que es un caso desbalanceado. Concluya sobre el resultado de la tabla del ANOVA obtenida. Nota: Hay tres formas fundamentalmente diferentes de ejecutar un ANOVA, se conocen como sumas de cuadrados de Tipo I, Tipo II y Tipo III. Los tres métodos dan el mismo resultado cuando el diseño está equilibrado. Sin embargo, cuando el diseño está desequilibrado, no dan los mismos resultados.
Planteamiento de hipotesis \[H_o: \mu_1 = \mu_2 = \mu_3 = \mu_4 = \mu_5\]
\[H_a: H_o \ es \ falsa\]
library(car)## Loading required package: carDataanova_1 <- aov(Pérdida ~ Granja * Método, data = df2)
Anova(anova_1, type = "III") #ANOVA para modelo desbalanceado con interaccion## Anova Table (Type III tests)
##
## Response: Pérdida
## Sum Sq Df F value Pr(>F)
## (Intercept) 55.311 1 6.3651 0.02071 *
## Granja 0.289 1 0.0333 0.85723
## Método 10.551 4 0.3035 0.87196
## Granja:Método 0.699 4 0.0201 0.99913
## Residuals 165.107 19
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Aunque no hay hipotesis planateadas con respecto a la interacción, al visualizar el efecto de los dos factores, se evidencia un p.valor de 0.99913, lo cual es un indicativo de que no hay presencia de interaccion.
Por lo tanto, se ejecutó el modelo de la siguiente forma / sin interacción:
library(car)
anova_2 <- aov(Pérdida ~ Granja + Método, data = df2)
Anova(anova_2, type = "III") #ANOVA para un diseño desbalanceado sin interaccion## Anova Table (Type III tests)
##
## Response: Pérdida
## Sum Sq Df F value Pr(>F)
## (Intercept) 222.230 1 30.8269 1.198e-05 ***
## Granja 0.085 1 0.0118 0.9146
## Método 47.794 4 1.6575 0.1941
## Residuals 165.806 23
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Apartir del analisis de varianza se puede concluir que el factor granja no genera mayor problema, debido a que el p.valor es mayor al 5%. Lo que sugiere que pudo haberse bloqueado o no, ya que esadisticamente no hubiese causado gran diferencia. Sin embargo, en el caso de los metodos este pvalor nos está indicando que estos no causan efectos significativos en la variable respuesta.
De manera que se comenzó con un diseño factorial simple en bloques con arreglo completamentente al azar con interacción, se pasó a un modelo sin interacción, y al ver que los bloques estadisticamente no generar mayor problema, se ejecutó un analisis de varianza para un diseño factorial simple:
library(car)
anova_4 <- aov(Pérdida ~ Método, data = df2)
Anova(anova_4, type = "III") #ANOVA para diseño factorial simple desbalanceado## Anova Table (Type III tests)
##
## Response: Pérdida
## Sum Sq Df F value Pr(>F)
## (Intercept) 432.26 1 62.5371 3.87e-08 ***
## Método 47.76 4 1.7275 0.1768
## Residuals 165.89 24
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1-¿Afecta el orden de colocación de los efectos del modelo dentro del software R? Verifique si la tabla del ANOVA cambia?.
library(car)
anova_3 <- aov(Pérdida ~ Método * Granja, data = df2)
Anova(anova_3, type = "III")## Anova Table (Type III tests)
##
## Response: Pérdida
## Sum Sq Df F value Pr(>F)
## (Intercept) 55.311 1 6.3651 0.02071 *
## Método 10.551 4 0.3035 0.87196
## Granja 0.289 1 0.0333 0.85723
## Método:Granja 0.699 4 0.0201 0.99913
## Residuals 165.107 19
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Modelo sin interacción:
library(car)
anova_5 <- aov(Pérdida ~ Método + Granja, data = df2)
Anova(anova_5, type = "III")## Anova Table (Type III tests)
##
## Response: Pérdida
## Sum Sq Df F value Pr(>F)
## (Intercept) 222.230 1 30.8269 1.198e-05 ***
## Método 47.794 4 1.6575 0.1941
## Granja 0.085 1 0.0118 0.9146
## Residuals 165.806 23
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1RTA:Al realizar la tabla de analisis de varianza de las dos maneras, se puede evidenciar que el orden de colocación de los efectos del modelo no afecta. Ya que, por esta razon se uso la suma de cuadrados tipo III,donde el ajuste incluye interacciones en niveles superiores. Porque, al calcular sumas de cuadrados para un diseño no balanceado de la misma manera que lo hace para un diseño balanceado generaría problemas. Debido a que a diferencia de un diseño equilibrado, las sumas de cuadrados para los efectos principales variarán según el orden en que se ingresen los factores al modelo. Esto se debe a que el diseño ya no es balanceado y, como resultado, algunas de las variables explicativas pueden estar correlacionadas. Por lo tanto, para los diseños no balanceados, es común utilizar sumas de cuadrados ajustadas que no se ven afectadas por el orden de los factores en el modelo (Tipo III).
B.Estimar el valor de la observación usando el promedio de los datos para los seis 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.
library(readxl) # Importar datos de excel
cotton <- read_excel("cotton.xlsx")data_8 <- data.frame(cotton) # Generar un data.frame
head(data_8)## Granja Metodo Perdida
## 1 1 M1 9.30
## 2 1 M2 5.54
## 3 1 M3 7.67
## 4 1 M4 7.89
## 5 1 M5 9.27
## 6 2 M1 6.75cotton$Metodo <- factor(cotton$Metodo)
cotton$Metodo## [1] M1 M2 M3 M4 M5 M1 M2 M3 M4 M5 M1 M2 M3 M4 M5 M1 M2 M3 M4 M5 M1 M2 M3 M4 M5
## [26] M1 M2 M3 M4 M5
## Levels: M1 M2 M3 M4 M5cotton$Granja <- factor(cotton$Granja)
cotton$Granja## [1] 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 4 4 4 4 4 5 5 5 5 5 6 6 6 6 6
## Levels: 1 2 3 4 5 6Planteamiento de hipotesis
\[H_o: \mu_1 = \mu_2 = \mu_3 = \mu_4 = \mu_5\]
\[H_a: H_o \ es \ falsa\]
• Factor Metodo:
Para evaluar el efecto de los tratamientos, la suma de cuadrados de tratamientos debe ajustarse por bloques, por lo tanto primero se introducen los bloques y después los tratamientos.
mod1 <- aov(Perdida~ Granja * Metodo, data = data_8)
summary(mod1)## Df Sum Sq Mean Sq F value Pr(>F)
## Granja 1 0.08 0.078 0.009 0.924
## Metodo 4 51.16 12.789 1.548 0.227
## Granja:Metodo 4 0.58 0.146 0.018 0.999
## Residuals 20 165.23 8.261Al visualizar el efecto de los dos factores, se puede evidenciar un p.valor de 0.999, lo cual nos permite afirmar que no hay presencia de interaccion.
Por lo tanto, se ejecutó el modelo de la siguiente forma:
mod2 <- aov(Perdida~ Granja + Metodo, data = data_8)
summary(mod2)## Df Sum Sq Mean Sq F value Pr(>F)
## Granja 1 0.08 0.078 0.011 0.916
## Metodo 4 51.16 12.789 1.851 0.152
## Residuals 24 165.81 6.909El p.valor 0.916 de la granja, sugiere que pudo haberse bloqueado o no, ya que esadisticamente no hubiese causado gran diferencia. Al igual que en el caso de los metodos.
Comparación con el diseño desbalanceado
library(car)
anova_2 <- aov(Pérdida ~ Granja + Método, data = df2)
Anova(anova_2, type = "III")## Anova Table (Type III tests)
##
## Response: Pérdida
## Sum Sq Df F value Pr(>F)
## (Intercept) 222.230 1 30.8269 1.198e-05 ***
## Granja 0.085 1 0.0118 0.9146
## Método 47.794 4 1.6575 0.1941
## Residuals 165.806 23
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod2 <- aov(Perdida~ Granja + Metodo, data = data_8)
summary(mod2)## Df Sum Sq Mean Sq F value Pr(>F)
## Granja 1 0.08 0.078 0.011 0.916
## Metodo 4 51.16 12.789 1.851 0.152
## Residuals 24 165.81 6.909Estadísticamente se hubiese llegado a las mismas conclusiones, ya que en ambos casos, tanto en el diseño balanceado como en el desbalanceado, el factor bloque pude no incluirse, al no ser significativo, ya que en ambos casos el p.valor es mayor 0.05, de igual manera ninguno de los anova proporciona sufiente informacion para tomar desciones sobre los metodos, porque estadísticamente tampoco son significativos. Por lo tanto se procedió a hacer un gráfico, del cual se puede evidenciar que el metodo en el que menos se pierde cantidad de algodón sería con el numero dos, por lo cual es el más recomendable.
library(ggplot2)
ggplot(data_8, aes(x = Metodo, y = Perdida, colour = Metodo))+
geom_boxplot() +theme_bw()
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.3+0.1 y de la capa inferior osciló entre 2 y 2.0+0.2. Use la función sort.int de R para ordenar los datos de cada capa con la opción partial=25+0 dentro de la propia función sort.int.
set.seed(2021) # Permite generar numeros aleatoreos
CO_5 = runif(50, min = 3.0, max = 3.4)
CO5 =sort.int(CO_5,partial = 25)
CO_10 = runif(50, min = 2, max = 2.2) # runif, permite generar datos de la distribucion uniforme
CO10 =sort.int(CO_10,partial = 25)
profundidad = gl(2, 50, 100, c('CO_5','CO_10')); profundidad # gl, usado para factores## [1] CO_5 CO_5 CO_5 CO_5 CO_5 CO_5 CO_5 CO_5 CO_5 CO_5 CO_5 CO_5
## [13] CO_5 CO_5 CO_5 CO_5 CO_5 CO_5 CO_5 CO_5 CO_5 CO_5 CO_5 CO_5
## [25] CO_5 CO_5 CO_5 CO_5 CO_5 CO_5 CO_5 CO_5 CO_5 CO_5 CO_5 CO_5
## [37] CO_5 CO_5 CO_5 CO_5 CO_5 CO_5 CO_5 CO_5 CO_5 CO_5 CO_5 CO_5
## [49] CO_5 CO_5 CO_10 CO_10 CO_10 CO_10 CO_10 CO_10 CO_10 CO_10 CO_10 CO_10
## [61] CO_10 CO_10 CO_10 CO_10 CO_10 CO_10 CO_10 CO_10 CO_10 CO_10 CO_10 CO_10
## [73] CO_10 CO_10 CO_10 CO_10 CO_10 CO_10 CO_10 CO_10 CO_10 CO_10 CO_10 CO_10
## [85] CO_10 CO_10 CO_10 CO_10 CO_10 CO_10 CO_10 CO_10 CO_10 CO_10 CO_10 CO_10
## [97] CO_10 CO_10 CO_10 CO_10
## Levels: CO_5 CO_10CO = c(CO_5,CO_10)
dataF = data.frame(CO, profundidad);
head(dataF)## CO profundidad
## 1 3.180507 CO_5
## 2 3.313512 CO_5
## 3 3.283873 CO_5
## 4 3.152698 CO_5
## 5 3.254530 CO_5
## 6 3.280538 CO_5AUse 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. ´
Ventanilla <- expand.grid(Longitud= seq(0,100,25), Latitud = seq(0,200,length.out = 10))
Ventanilla## Longitud Latitud
## 1 0 0.00000
## 2 25 0.00000
## 3 50 0.00000
## 4 75 0.00000
## 5 100 0.00000
## 6 0 22.22222
## 7 25 22.22222
## 8 50 22.22222
## 9 75 22.22222
## 10 100 22.22222
## 11 0 44.44444
## 12 25 44.44444
## 13 50 44.44444
## 14 75 44.44444
## 15 100 44.44444
## 16 0 66.66667
## 17 25 66.66667
## 18 50 66.66667
## 19 75 66.66667
## 20 100 66.66667
## 21 0 88.88889
## 22 25 88.88889
## 23 50 88.88889
## 24 75 88.88889
## 25 100 88.88889
## 26 0 111.11111
## 27 25 111.11111
## 28 50 111.11111
## 29 75 111.11111
## 30 100 111.11111
## 31 0 133.33333
## 32 25 133.33333
## 33 50 133.33333
## 34 75 133.33333
## 35 100 133.33333
## 36 0 155.55556
## 37 25 155.55556
## 38 50 155.55556
## 39 75 155.55556
## 40 100 155.55556
## 41 0 177.77778
## 42 25 177.77778
## 43 50 177.77778
## 44 75 177.77778
## 45 100 177.77778
## 46 0 200.00000
## 47 25 200.00000
## 48 50 200.00000
## 49 75 200.00000
## 50 100 200.00000profundidad <- seq(0, 100, l = 100)
CO <- seq(0, 200, l = 200)
pred.grid <- expand.grid(x = profundidad, y = CO)
plot(dataF$CO, pch = 20)
points(pred.grid, pch = 3, cex = 0.2)
B.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.
dataor <- data.frame(Longitud = rep(Ventanilla$Longitud,2),Latitud = rep(Ventanilla$Latitud,2),profundidad = rep(c(5,10),each = 50), carbono = c(CO_5,CO_10))
library(plotly)##
## 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':
##
## layoutp <- plot_ly(dataor, x = dataor$Longitud, y = dataor$Latitud, z = dataor$profundidad, color = dataor$carbono, colors = c("#ffce99",'#ccaf8f', "#e09867", '#6e5538')) %>% add_markers() %>% layout( scene = list(xaxis = list(title = 'Longitud'), yaxis = list(title = 'Latitud'), zaxis = list(title = 'Profundidad'))); p## 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.library(scatterplot3d)
scatterplot3d( dataF$profundidad, dataF$CO, xlab = "l", ylab = "f", zlab = "profundidad") 
D.Compare si se encuentran diferencias en la media de carbono entre capas utilizando un nivel de confianza del 95%.
Planteamiento de Hipotesis:
\[H_o:\mu_{CO_5} = \mu_{CO_{10}} \] \[H_a:\mu_{CO_5} \neq \mu_{CO_{10}} \]
Para esto se hace uso de la prueba t-para dos muestras dependientes/pareadas.
Aplicando la prueba t para muestras pareadas:
pruebat = t.test(CO_5, CO_10, alternative ="t",mu = 0,paired = TRUE,conf.level = 0.95)
pruebat##
## Paired t-test
##
## data: CO_5 and CO_10
## t = 67.658, df = 49, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 1.110511 1.178499
## sample estimates:
## mean of the differences
## 1.144505ifelse(pruebat$p.value<0.05, 'Rechazo Ho', 'No rechazo Ho')## [1] "Rechazo Ho"Por lo tanto se rechaza la hipotesis nula utilizando un nivel de confianza del 95%, la cual dicta que no se encuentran diferencias en la media de carbono entre capas, por lo tanto los datos generados no se comportan estadísticamente iguales, y hay diferencias de los niveles de carbono con la profundidad.
4.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:
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 repeticones por tratamiento
set.seed(2020)
D <- D[order(sample(1:18)),]
class(D)## [1] "data.frame"D$biomasa= sort.int(rnorm(18,3,0.3),partial = 9) # Crea la respuesta
head(D)## F1 F2 biomasa
## 2 3.75 4 2.708826
## 10 3.25 4 2.772692
## 16 3.25 6 2.143359
## 4 3.25 5 2.560519
## 13 3.25 5 2.708666
## 6 4.25 5 2.773705A.Escriba (completamente especificado) el modelo del diseño
\[y_{ijk} = \mu + \tau_i + \delta_j + (\tau\delta)_{ij}+ \epsilon_{ijk} \\ i: 1..3 \\ j: 1..2 \\ k: 1..2\\ Ademas~de~condiciones~laterales~respectivas\] Donde:
\(\mu\), es la media general
\(\tau_i\), es el efecto de \(i-esimo~nivel~del~factor~A\)
\(\delta_i\), es el efecto de \(j-esimo~nivel~del~factor~B\)
\((\tau\delta)_{ij}\), efecto causado por la interacción del nivel \(i\) del factor A y el nivel \(j\) del factor B.
\(\epsilon_{ijk}\), termino de error
B.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.Además , use los resultados del ANOVA y el gráfico de interacción (ggplot2) para visualizar si existe o no interacción entre los factores.
dosis = factor(D$F1, labels = c("3.25", "3.75", "4.25"))
aplicacion = factor(D$F2, labels = c(4, 5, 6))Para realizar el primer análisis de varianza se toman en cuenta los factores que en este caso son la dosis del insecticida y el número de aplicaciones en el desarrollo del cultivo y la posible interacción que pueda existir entre entre ellos.
mod5 = aov(D$biomasa~dosis*aplicacion)
summary(mod5)## Df Sum Sq Mean Sq F value Pr(>F)
## dosis 2 0.4161 0.20804 1.979 0.194
## aplicacion 2 0.3426 0.17128 1.630 0.249
## dosis:aplicacion 4 0.1635 0.04087 0.389 0.812
## Residuals 9 0.9459 0.10510Al observar los p.valores, se puede ver que los factores no tienen efectos significativos.También, se puede ver que la interacción tampoco, por tanto, se realiza un análisis de varianza , descartando la interacción entre los factores.
mod5 = aov(D$biomasa~dosis+aplicacion)
summary(mod5)## Df Sum Sq Mean Sq F value Pr(>F)
## dosis 2 0.4161 0.20804 2.438 0.126
## aplicacion 2 0.3426 0.17128 2.007 0.174
## Residuals 13 1.1094 0.08534De nuevo, se puede observar en la tabla del análisis de varianza que los efectos de los factores no están siendo significativos.Sin embargo se deben revisar los supuestos para poder validar los análisis que se han hecho basados en el ANOVA
1.Normalidad de los residuos del modelo
resid5 = mod5$residuals
norm5 = shapiro.test(resid5)
ifelse(0.05>norm5$p.value , "NO_normales" , "Residuos_normales")## [1] "Residuos_normales"- Igualdad de varianzas entre tratamientos
bar = bartlett.test(resid5 ~ interaction(aplicacion,dosis))
ifelse(0.05>bar$p.value , "Var_Desiguales" , "Varianzas_iguales")## [1] "Varianzas_iguales"- Independencia de los residuos
plot(resid5)
En el grafico se evidencia independencia entre los residuos, ya que no siguen un patrón en especifíco.
C El investigador quiso colocar como covariable el contenido de arcilla(expansible) en el suelo utilizado en cada unidad experimental. Use los resultados del ANOVA y el gráfico de interacción (ggplot2) para visualizar si existe o no interacción entre los factores.
Dosis = factor(D$F1, labels = c("3.25", "3.75", "4.25"))
Aplicacion = factor(D$F2, labels = c(4, 5, 6))
datao = data.frame(Dosis,Aplicacion)
datao$biomasa = D$biomasa
head(datao)## Dosis Aplicacion biomasa
## 1 3.75 4 2.708826
## 2 3.25 4 2.772692
## 3 3.25 6 2.143359
## 4 3.25 5 2.560519
## 5 3.25 5 2.708666
## 6 4.25 5 2.773705library(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:car':
##
## recode## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionattach(datao)## The following objects are masked _by_ .GlobalEnv:
##
## Aplicacion, Dosisz<-aggregate(x=list(media_B=datao$biomasa), by=list(Dosis=datao$Dosis,Aplicacion=datao$Aplicacion),
FUN=mean, na.rm=TRUE)
z## Dosis Aplicacion media_B
## 1 3.25 4 2.965294
## 2 3.75 4 3.080186
## 3 4.25 4 3.344110
## 4 3.25 5 2.634592
## 5 3.75 5 3.128949
## 6 4.25 5 2.894908
## 7 3.25 6 2.592757
## 8 3.75 6 2.801410
## 9 4.25 6 3.021714library(ggplot2)
z %>%
ggplot(aes(Dosis,media_B,color=Aplicacion))+
geom_line(aes(group = Aplicacion))
interaction.plot(datao$Dosis,datao$Aplicacion,datao$biomasa)
En los dos gráficos de interacción se puede envidenciar una posible presencia de este fenómeno entre los dos factores.Sin embargo , el ANOVA cuyos supuestos fueron evaluados, nos permite decir que esta interacción no es estadísticamente significativa.Además de estos gráficos de interacción podemos concluir que Estadísticamente el mejor tratamiento es de una dosis de 4.25 con 4 aplicaciones.Sin embargo , desde el punto de vista agronómico seria recomendable evaluar los costos y los impactos que tendría cada tratamiento.
D.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.
set.seed(2020)
arcilla = sort.int(runif(18, 0.20, 0.40)); arcilla## [1] 0.2005165 0.2134769 0.2258305 0.2272194 0.2786236 0.2788452 0.2818575
## [8] 0.2845458 0.2953782 0.3079385 0.3237004 0.3240412 0.3293806 0.3307115
## [15] 0.3487672 0.3528828 0.3652331 0.3921445E = data.frame (D, arcilla);
head(E)## F1 F2 biomasa arcilla
## 2 3.75 4 2.708826 0.2005165
## 10 3.25 4 2.772692 0.2134769
## 16 3.25 6 2.143359 0.2258305
## 4 3.25 5 2.560519 0.2272194
## 13 3.25 5 2.708666 0.2786236
## 6 4.25 5 2.773705 0.2788452dosiss = factor(E$F1, labels = c("3.25", "3.75", "4.25"))
aplicaciones = factor(E$F2, labels = c(4, 5, 6))
ancova1 = aov(biomasa~ arcilla+dosiss*aplicaciones,data=E )
summary(ancova1)## Df Sum Sq Mean Sq F value Pr(>F)
## arcilla 1 0.8839 0.8839 21.347 0.00171 **
## dosiss 2 0.1496 0.0748 1.806 0.22526
## aplicaciones 2 0.4777 0.2388 5.768 0.02812 *
## dosiss:aplicaciones 4 0.0256 0.0064 0.155 0.95541
## Residuals 8 0.3313 0.0414
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Como la interacción no fue significativa en el ancova que se realizó, se procede a correr un modelo aditivo (sin interacción)
ancova2 = aov(biomasa~ arcilla+(dosiss+aplicaciones),data=E )
summary(ancova2)## Df Sum Sq Mean Sq F value Pr(>F)
## arcilla 1 0.8839 0.8839 29.721 0.000147 ***
## dosiss 2 0.1496 0.0748 2.515 0.122440
## aplicaciones 2 0.4777 0.2388 8.030 0.006117 **
## Residuals 12 0.3569 0.0297
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Con los resultados de este ANCOVA , se puede decir que el factor aplicaciones y la covariable está teniendo efectos significagtivos en la variable respuesta (Biomasa). Sin embargo, se deben probar los supuestos de este modelo para poder evaluar los efectos individuales del factor “aplicaciones”.
Probando los supuestos del ANCOVA
1.Normalidad de los residuos del modelo
residanc = ancova2$residuals
normanc = shapiro.test(residanc)
ifelse(0.05>normanc$p.value , "NO_normales" , "Residuos_normales")## [1] "Residuos_normales"2. Igualdad de varianzas entre tratamientos
baranc = bartlett.test(residanc ~ interaction(aplicaciones,dosiss))
ifelse(0.05>bar$p.value , "Var_Desiguales" , "Varianzas_iguales")## [1] "Varianzas_iguales"3. Independencia de los residuos
plot(residanc)
4.Existe una relación lineal entre la covariable (Arcilla) y la variable respuesta (Biomasa)
cor.test(E$arcilla, E$biomasa)##
## Pearson's product-moment correlation
##
## data: E$arcilla and E$biomasa
## t = 3.7909, df = 16, p-value = 0.001603
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.3255678 0.8740494
## sample estimates:
## cor
## 0.6878814Este supuesto se estaría cumpliendo, pues en el test de correlación se puede ver un p valor menor al 5%
5. Homogeneidad de las pendientes de regresión
library(rstatix)##
## Attaching package: 'rstatix'## The following object is masked from 'package:stats':
##
## filterE %>% anova_test(biomasa ~ arcilla + F1 + F2 + F1*F2 + arcilla*F1 + arcilla*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 19.052 0.001 * 0.656
## 2 F1 1 10 0.167 0.691 0.016
## 3 F2 1 10 13.148 0.005 * 0.568
## 4 F1:F2 1 10 0.165 0.693 0.016
## 5 arcilla:F1 1 10 3.437 0.093 0.256
## 6 arcilla:F2 1 10 0.033 0.859 0.003
## 7 arcilla:F1:F2 1 10 0.863 0.375 0.079Para probar este supuesto, se realizo una ultima prueba donde se pudiera observar la interacción de la covariable con los factores , de este analisis se puede concluir que hubo homogeneidad de las pendientes de regresión, ya que los términos de interacción entre las variables de covariable (arcilla) y de los factores (dosis y aplicaciones) no fueron estadísticamente significativos, p > 0,05.
Al cumplirse todos los supuestos del ANCOVA realizado se puede confirmar que la presencia de arcillas infuye en la respuesta de la biomasa.Aunque , según el coeficiente de correlación de pearson encontrado en el test de correlación se puede concluir que estas no afectan de manera negativa a la variable respuesta que en este caso es biomasa, pues el valor de este es positivo.Además, se pueden evaluar los efectos indivivuales de un factor a través de un test Tukey de comparaciones múltiples.
t1 = TukeyHSD(ancova2, 'aplicaciones');t1## Warning in replications(paste("~", xx), data = mf): non-factors ignored: arcilla## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = biomasa ~ arcilla + (dosiss + aplicaciones), data = E)
##
## $aplicaciones
## diff lwr upr p adj
## 5-4 -0.2248984 -0.4905335 0.04073678 0.1011482
## 6-4 -0.3936393 -0.6592745 -0.12800419 0.0050315
## 6-5 -0.1687410 -0.4343761 0.09689418 0.2469676Como se puede observar en los gráficos para intervalos de confianza en el test de TUkey y un p valor menor al 5% se puede decir que las comparaciones entre los niveles 6-4 de las aplicaciones son estaditícamente significativos.
plot(t1 , las=1 , col="brown")
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.
knitr::include_graphics("/cloud/project/3.png") 
Figura 1 Diseño anidado escalonado.
library(readxl)
escalonado <- read_excel("escalonado.xlsx")
datos_p = data.frame(escalonado)
head(datos_p)## Finca Variedad Test Parcela Respuesta
## 1 A 1 1 1 9.76
## 2 A 1 1 2 10.65
## 3 A 1 1 3 6.50
## 4 A 1 1 4 8.08
## 5 A 1 1 5 1.84
## 6 A 1 1 6 9.00A.Use la librería lme4 tal como aparece en el código que aparece 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.
library(daewr)## Registered S3 method overwritten by 'DoE.base':
## method from
## factorize.factor conf.design anova_p<-aov(Respuesta ~ Parcela + Parcela:Finca + Parcela:Finca:Variedad, data = datos_p)
summary(anova_p)## Df Sum Sq Mean Sq F value Pr(>F)
## Parcela 1 2.8 2.755 0.299 0.586
## Parcela:Finca 1 3.7 3.674 0.398 0.530
## Parcela:Finca:Variedad 1 6.3 6.279 0.681 0.412
## Residuals 76 700.9 9.222Se observa que ninguno tiene efecto significativo, debido a que todos los p.valores son mayores al 0.05, por lo tanto el analisis de varianza no aporta infromación para poder tomar una determinación con respecto al experimento, de manera que se realizó un analisis descriptivo, por medio de un boxplot, para ver como se comportar la variedad dentro de cada finca.
ggplot(data = datos_p, aes(x = Variedad, y = Respuesta, colour = Finca,)) +
geom_boxplot() +
theme_bw()
B.Use los datos que se muestran para estimar las varianzas antes descritas. Una ayuda para la solución de este problema puede encontrarse en el libro: Design and Analysis of Experiments with R de John Lawson.
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':
##
## cakemodr3 <- lmer( Respuesta ~ 1 + (1|Parcela) + (1|Parcela:Finca)+
+ (1|Parcela:Finca:Variedad), data = datos_p)## boundary (singular) fit: see ?isSingularsummary(modr3)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## Respuesta ~ 1 + (1 | Parcela) + (1 | Parcela:Finca) + +(1 | Parcela:Finca:Variedad)
## Data: datos_p
##
## REML criterion at convergence: 335.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.48350 -0.37362 0.01905 0.46662 1.75836
##
## Random effects:
## Groups Name Variance Std.Dev.
## Parcela:Finca:Variedad (Intercept) 1.166 1.080
## Parcela:Finca (Intercept) 0.000 0.000
## Parcela (Intercept) 7.167 2.677
## Residual 1.176 1.085
## Number of obs: 80, groups:
## Parcela:Finca:Variedad, 60; Parcela:Finca, 40; Parcela, 20
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 8.1680 0.6273 13.02
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular6.En el enlace https://cran.r-project.org/web/packages/asbio/asbio.pdf se tienen unos datos de potasio de muestras de suelos medidas en 8 diferentes laboratorios.
library(readxl)
datak <- read_excel("datak.xlsx")
datak = data.frame(datak)
lab = factor(datak$lab); lab## [1] B B B B B B B B B D D D D D D D D D E E E E E E E E E F F F F F F F F F G G
## [39] G G G G G G G H H H H H H H H H I I I I I I I I I J J J J J J J J J
## Levels: B D E F G H I Jdatos = data.frame(lab, soil = datak$K);
head(datos)## lab soil
## 1 B 296
## 2 B 260
## 3 B 341
## 4 B 359
## 5 B 323
## 6 B 321A.Compare descriptivamente (medidas, tablas y gráficos) para representar los datos. #### Estadística descriptiva.
library(psych)##
## Attaching package: 'psych'## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha## The following object is masked from 'package:car':
##
## logitdescribe(datos$soil, ranges=T)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 72 307.79 48.4 314 311.1 43.74 187 413 226 -0.54 0.04 5.7describeBy(datos$soil, datos$lab)##
## Descriptive statistics by group
## group: B
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 9 326.11 44.41 323 326.11 40.03 260 413 153 0.41 -0.7 14.8
## ------------------------------------------------------------
## group: D
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 9 321.11 29.26 326 321.11 25.2 266 354 88 -0.68 -1.02 9.75
## ------------------------------------------------------------
## group: E
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 9 316.56 80.35 351 316.56 72.65 187 400 213 -0.54 -1.44 26.78
## ------------------------------------------------------------
## group: F
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 9 315.67 25.05 308 315.67 23.72 274 354 80 0.05 -1.22 8.35
## ------------------------------------------------------------
## group: G
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 9 304.11 15.37 301 304.11 22.24 280 326 46 -0.14 -1.51 5.12
## ------------------------------------------------------------
## group: H
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 9 229.33 23.89 226 229.33 22.24 199 280 81 0.74 -0.32 7.96
## ------------------------------------------------------------
## group: I
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 9 313.22 31.4 311 313.22 41.51 269 355 86 -0.09 -1.81 10.47
## ------------------------------------------------------------
## group: J
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 9 336.22 21.61 342 336.22 20.76 299 359 60 -0.46 -1.53 7.2mean<-tapply(datos$soil, datos$lab,mean);mean## B D E F G H I J
## 326.1111 321.1111 316.5556 315.6667 304.1111 229.3333 313.2222 336.2222sd<-tapply(datos$soil, datos$lab,sd);sd## B D E F G H I J
## 44.40564 29.25510 80.35097 25.04995 15.37404 23.89037 31.39577 21.60890median<- tapply(datos$soil, datos$lab,median); median## B D E F G H I J
## 323 326 351 308 301 226 311 342library(ggplot2)
ggplot(datos, aes(x = lab, y = soil,colour = lab))+
geom_boxplot()
B.¿qué prueba me recomendaría para comparar la medida que usted seleccione? Proponga una solución. Sabiendo que son muestras mezcladas de una misma finca, ¿ Se perciben diferencias en las medidas como consecuencia probable de los laboratorios? Sugerencia: Use el enfoque no paramétrico considerado en clase y su respectiva prueba de comparación por pares (Nemenyi)
Como se observa a continuación, en primera instancia se realizo el anova, pero al verificar los supuestos se evidenció de que no se cumple la normalidad de los residuales ni la homosedasticidad de las varianzas, por lo cual se recomienda realizar el test de Kruskal-Wallis.
mod2<-aov(datos$soil~datos$lab)
summary(mod2)## Df Sum Sq Mean Sq F value Pr(>F)
## datos$lab 7 68930 9847 6.472 8.92e-06 ***
## Residuals 64 97370 1521
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(mod2, 1)
plot(mod2, 2)
Prueba de normalidad del vector de resiudales del modelo:
\[H_0:Los~residuales~poseen~distribucion~normal\] \[H_0:Los~residuales~no~poseen~distribucion~normal\]
resid = mod2$residuals
shapiro.test(resid)##
## Shapiro-Wilk normality test
##
## data: resid
## W = 0.94688, p-value = 0.004268# Claramente se observa que los residuales no son normales (p-value < 0.05)Prueba de igualdad de varianza:
\[H_0: \sigma^2_{{B}}=\sigma^2_{{C}}=\sigma^2_{{D}}=\sigma^2_{{F}}=\sigma^2_{{G}}=\sigma^2_{{H}}=\sigma^2_{{I}}=\sigma^2_{{J}}\]
bartlett.test(resid~datos$lab)##
## Bartlett test of homogeneity of variances
##
## data: resid by datos$lab
## Bartlett's K-squared = 32.201, df = 7, p-value = 3.727e-05# Claramente no se observa homosedasticidad (p-value < 0.05)Test de Kruskal-Wallis
\[H_0:\mu_{rangoB}=\mu_{rangoD}=\mu_{rangoF}=\mu_{rangoG}=\mu_{rangoH}=\mu_{rangoI}=\mu_{rangoJ}\]
kruskal.test(datos$soil~datos$lab) ##
## Kruskal-Wallis rank sum test
##
## data: datos$soil by datos$lab
## Kruskal-Wallis chi-squared = 24.482, df = 7, p-value = 0.000937qchisq(0.05, 7, lower.tail = F)## [1] 14.06714Al realizar la comparacion, se evidencia que el test encuentra diferencias de al menos dos grupos, es decir que los rangos promedios asignados a cada grupo son diferentes. De manera que, esta situación podria deberse a causa de que las muestras fueron tratadas en distintos laboratorio. Por ello, se ejecutó la prueba de comparación por pares (Nemenyi) para corroborar:
PMCMR::posthoc.kruskal.nemenyi.test(datos$soil~datos$lab)##
## Pairwise comparisons using Tukey and Kramer (Nemenyi) test
## with Tukey-Dist approximation for independent samples
##
## data: datos$soil by datos$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: noneDesde una manera general, se puede envidenciar que se encuentran diferencias entre los resultados del algunos laboratorios, pero evidentemente no es el caso de todos.
Independencia de los residuales
plot(mod2$residuals)
7.Diseñe un experimento en parcelas divididas en bloques completos (diseño en franjas o strip plot design). Genere los datos usted mismo y esquematice el diseño.
Hay 6 campos enteros, el primer factor que se establece de forma vertical corresponde a tres niveles de variedades de arroz (V1, V2 y V3).Posteriomente,se trazan franjas horizontales donde se le asignó aleatoriamente uno de los tres niveles de fertilización. Se determinó el rendimiento de arroz (variable respuesta).
strip <- read.csv2("/cloud/project/strip.csv")
head(strip)## BLOQUE Variable Fertilizante yield
## 1 1 V1 0 103.60
## 2 1 V1 60 126.25
## 3 1 V1 120 115.39
## 4 1 V2 0 103.80
## 5 1 V2 60 119.32
## 6 1 V2 120 106.40library(collapsibleTree)
collapsibleTree(strip,hierarchy = c("BLOQUE", "Variable", "Fertilizante"),zoomable = T, fill = "#ceebe1")A.Expliqué las razones de colocar el primer factor en la parcela principal y el segundo en la subparcela.
En primer lugar, se debe resaltar que el orden depende del investigador. En el caso de las variedades de arroz, se tuvo en cuenta que el factor más dificil de aleatorizar en campo, razon por la cual los niveles de fertilizante fueron aleatorizadas entre las variedades. Por lo tanto, la parcela principal vendría siendo la variedad, y la subparcela sería la fertilización. Por cuestiones netamente experimentales y de campo.
strip$Variable = as.factor(strip$Variable)
strip$Fertilizante = as.factor(strip$Fertilizante)attach(strip)
library(agricolae)
anovastrip = strip.plot(BLOCK = BLOQUE, COL = Variable, ROW = Fertilizante, Y = yield)##
## ANALYSIS STRIP PLOT: yield
## Class level information
##
## Variable : V1 V2 V3
## Fertilizante : 0 60 120
## BLOQUE : 1 2 3 4 5 6
##
## Number of observations: 54
##
## model Y: yield ~ BLOQUE + Variable + Ea + Fertilizante + Eb + Fertilizante:Variable + Ec
##
## Analysis of Variance Table
##
## Response: yield
## Df Sum Sq Mean Sq F value Pr(>F)
## BLOQUE 5 1245.73 249.15 4.0425 0.01068 *
## Variable 2 868.93 434.47 1.9446 0.19348
## Ea 10 2234.27 223.43 3.6251 0.00686 **
## Fertilizante 2 427.02 213.51 1.9988 0.18609
## Eb 10 1068.19 106.82 1.7332 0.14166
## Fertilizante:Variable 4 219.83 54.96 0.8917 0.48704
## Ec 20 1232.65 61.63
## ---
## 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.1231Covariable
Como covariable se estableció la cantidad de un mineral , que se cree influye de manera negativa en el rendimiento del cultivo de arroz
set.seed(2021)
Cova = runif(54, min = 53 , max = 95); Cova## [1] 71.95323 85.91875 82.80665 69.03326 79.72560 82.45653 79.89843 64.20055
## [9] 87.24770 94.28545 54.14522 88.17460 78.33611 76.83304 87.44222 63.56599
## [17] 74.23075 89.43660 93.24363 75.91930 58.86239 93.12425 69.48474 64.27638
## [25] 77.03269 91.31004 92.24027 89.98042 92.71918 87.22955 54.37659 92.59379
## [33] 92.80501 90.88783 76.19518 62.44522 92.08995 63.08275 92.26582 71.78802
## [41] 60.26528 88.55961 70.90653 93.20682 58.68769 86.34201 90.23294 53.70075
## [49] 55.18769 87.57983 75.76241 80.61310 56.05648 75.31045strip2 = data.frame (strip, Cova);
head(strip2)## BLOQUE Variable Fertilizante yield Cova
## 1 1 V1 0 103.60 71.95323
## 2 1 V1 60 126.25 85.91875
## 3 1 V1 120 115.39 82.80665
## 4 1 V2 0 103.80 69.03326
## 5 1 V2 60 119.32 79.72560
## 6 1 V2 120 106.40 82.45653Antes de realizar el análisis de covarianza , se quiere evaluar uno de los supuestos que debe cumplir la covariable (Tener una relación de forma lineal con la variable respuesta) mediante es test de pearson para la correlación
cor.test(Cova,yield)##
## Pearson's product-moment correlation
##
## data: Cova and yield
## t = 0.011384, df = 52, p-value = 0.991
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.2662944 0.2692255
## sample estimates:
## cor
## 0.001578743En la salida del código, se puede observar que la variable respuesta no tiene correlación con la covariable que se planteó inicialmente , por tanto , se incumpliria uno de los supuestos del ANCOVA y no e podrían concluir sobre las respuestas que generaría este análisis
8. Realice un resumen con la nota que aparece en las siguientes direcciones sobre:
A.El uso de los diseños en parcelas divididas: http://207.67.83.164/quality-progress/2007/10/laboratory/when-should-you-consider-a-split-plot-design.html
B.Sobre lo que significa unidad experimental y unidad de observación Short communication: On recognizing the proper experimental unit in animal studies in the dairy sciences -https://www.sciencedirect.com/science/article/pii/S002203021630621X -https://online.stat.psu.edu/stat502/lesson/6/6.1-0
C.Guía para diseñar experimentos exitosos https://acsess.onlinelibrary.wiley.com/doi/full/10.2134/agronj2013.0114
Resumen. Inicialmente, en el primer link, titulado “When Should You Consider A Split-Plot Design?”, realizado por Christine M. Anderson-Cook, se aborda el diseño de parcelas dividías, el cual debe considerarse cuando en el experimento existan algunos factores con niveles más difíciles de cambiar (factores de parcela completa) que otros (factores de subparcelas). En relación a los factores difíciles de cambiar, en muchos casos son más costosos o requieren mucho tiempo, entonces es viable seleccionar un diseño que contenga menos parcelas completas, lo que representaría un ahorro potencialmente significativo. Sin embargo, es importante considerar qué significa “óptimo” para un experimento y centrarse en los criterios más relevantes para los objetivos del proyecto. De hecho, se debe ser muy cuidadoso, al momento de analizar un diseño de parcelas divididas, ya que puede darse una confusión y abordarlo como un diseño completamente aleatorizado, lo cual no es correcto, y podría llevar a conclusiones erróneas en el experimento. Otro aspecto interesante, abordado por Anderson, es que, en la determinación del orden de ejecución del experimento de parcelas divididas, hay dos aleatorizaciones separadas, en la primera se determina el orden en el que ejecuta las parcelas completas; y luego se recopilan observaciones dentro de cada parcela completa. Lo cual, conduce a dos términos de error en el modelo, lo que constituye una ventaja de estos diseños, ya que los parámetros del modelo medio se pueden estimar por separado de los términos de error de la parcela completa y de la subparcela. Asimismo, se menciona, que cuando los niveles de factor no se restablecen para ejecuciones consecutivas que requieran el mismo nivel de un factor, se conoce como parcelas divididas invertidas.
En segunda instancia, en el artículo titulado “Short communication: On recognizing the proper experimental unit in animal studies in the dairy sciences”, se aborda principalmente aspectos relacionados con el diseño y análisis estadístico de experimentos en lácteos, junto con la preocupación relacionada a las políticas implementadas en las revistas de investigación animal, en relación a las publicaciones. Ya que, para poder tener datos científicos replicables y reproducibles, es primordial haber implementado adecuadamente la estadística en el análisis de los datos. Para empezar, en los primeros párrafos, se hace especial distinción entre unidad experimental (a la que se aplica un tratamiento de forma independiente) y unidad de observación (entidad física sobre la que se mide un resultado de interés en un experimento), haciendo alusión a que puede llegar a confundirse los términos, y por ende afectar inferencias y análisis en el proceso. De hecho, la especificación de la unidad experimental no es cuestión de opinión; está determinado por la forma como se organizó el experimento, cómo se recopilaron los datos y el alcance previsto de la inferencia. En relación a la acción de inferir, según el artículo es muy relevante hacerlo, debido a que conlleva a una aplicabilidad más amplia de las conclusiones. De manera que, el alcance adecuado de la inferencia debe ser reconocido, con el fin de proporcionar un contexto para interpretar los resultados de cualquier estudio, de forma que se puedan abordar las circunstancias generales y específicas para la reproducibilidad de los resultados de la investigación. Finalmente, un concepto que desconocía era el de estructura de datos jerárquica, y según el artículo se refiere a una configuración de los datos donde se realizan las observaciones, las cuales no son independientes entre sí, sino que tienen una estructura de correlación impuesta por el diseño experimental. De hecho, “What Would Fisher Do?” (WWFD), es una estrategia general que, en el contexto de datos asumidos como normales, es esencial para identificar el error experimental adecuado y así distinguir entre una unidad experimental y una unidad de observación.
Asimismo, en el apartado 6.1 de la Universidad de Pensilvania, se define a la unidad experimental como aquella que recibe el tratamiento. Para ilustrar esto, se proporciona un ejemplo, dónde se buscaba conocer el efecto del agua de un arroyo contaminado en las lesiones de los peces. Para ello, se instalaron 2 acuarios (uno con agua contaminada vs. control), cada uno con 50 peces. La confusión podría radicar en considerar a los peces como unidades experimentales, cuando en realidad los niveles de tratamiento de agua se aplicaron a todo el acuario, por lo que la unidad experimental es un acuario con 50 peces. Adicionalmente, se abordó el concepto, que constituye a la réplica de un experimento, lo que se refiere a cada vez que se aplica el conjunto completo de niveles de tratamientos. Para el caso del ejemplo no hay replica alguna, y esto se conoce como “estudio no replicado”. Finalmente, el no reconocer las unidades de muestreo puede resultar un gran problema, conocido como la pseudo-replicación, donde cada unidad de muestreo es tratada como si fuera una unidad experimental , lo que consecuentemente inflaría el Error df, reduciría el MSE y aumentaría el valor del estadístico F (incorrecto). Finalmente, en el texto “Fundamentals of Experimental Design: Guidelines for Designing Successful Experiments”, se abordan una serie de consejos para reducir la probabilidad de fracaso en los experimentos, ya que el fracaso tiene grandes implicaciones, tanto económicas como emocionales. Para ello, el autor se basa en los principios básicos del diseño experimental: replicación, aleatorización, bloqueo y tamaño de las unidades experimentales. Especialmente, se hace énfasis en los experimentos biológicos, donde el fracaso no es una opción. Ya que, en este ámbito, se realizan mayoritariamente experimentos comparativos, entre sistemas que pueden tener relevancia para la investigación y avances científicos. De hecho, contar con una base de datos producto de largos periodos de investigación, puede ser muy útil, ya que puede ser usada con fines de planificación, para derivar estimaciones de las variaciones y expectativas de las diferencias medias de tratamiento, lo que permite al individuo realizar análisis de poder significativos. Sim embargo, cuando no se tengan acceso a estas bases de datos, el autor recomienda, comenzar a construirla. Asimismo, el investigador debe realizar un buen diseño experimental, junto con un previo y posterior análisis de datos. Y en el caso de que sean experimentos de grandes proporciones, se sugiere usar diseños de bloques incompletos, los cuales sean lo suficientemente flexibles para generar bloques organizados. Finalmente, un consejo muy importante, es el de extender los límites de la imaginación, base de conocimientos y experiencias.
9.Seleccionar un artículo científico de una revista de agronomía donde se haya utilizado un diseño en parcelas divididas. Hacer las críticas constructivas sobre:
A.La mención de la estructura factorial:
Diseño experimental de bloques al azar en parcelas divididas con 3 repeticiones.
-Fecha de siembra como factor principal (Ambiente) -Cultivar como subfactor. Se usaron cinco cultivares: DM2800 - A3901 - DM4800 - A5409 - A5901
La variable respuesta, es el rendimiento del cultivo de soja, el cual fue medido como materia seca por hectárea. De hecho, la materia seca por hectárea se obtiene a partir de un proceso de secado de la biomasa obtenida.
B.La razón de colocar cada factor en la parcela principal o en la subparcela:
El cultivar es más sencillo de aleatorizar, ya que no habría una restricción temporal, como sí ocurre en el caso de las fechas de siembra
C.La revisión de supuestos para el análisis de varianza:
Para comprobar los supuestos del modelo Anova se llevaron a cabo diferentes tests. En primera instancia, para verificar la normalidad de los errores se utilizó el test de Shapiro-Wilk, donde no se rechazó la hipótesis nula de normalidad. Para la independencia de los residuos se usó el gráfico que se muestra acontinuación
Cuantiles normales vs residuales obtenidos del modelo
Además, se aplicó el test de Breusch-Pagan para heteroscedasticidad. El resultado muestra que se rechaza el supuesto de homocedasticidad.
Distribución de los Residuos vs Valores estimados del Modelo
Se evidencia que la heteroscedastidad de los residuales, es decir que los errores no tienen varianza constante, por lo tanto, en el articulo se realiza una corrección por heteroscedasticidad en un segundo modelo aplicado. Para ello se estimó un modelo semejante al Anova con la identificación de cluster o agrupamiento en cada repetición. Sin embargo, el supuesto evaluado en el primer ANOVA no se estaría cumpliendo, pero aun así los investigadores usan los resultados para compararlos con los modelos lineales, y concluyen apartir de ambas situaciones.
D La tabla del análisis de varianza:
Como se evidencia en la tabla anterior, hay presencia de interacción entre los dos factores en estudio (ambiente y cultivar), con un p.valor inferior a 0.05.
E El uso de muchos análisis de varianzas en lugar de uno solo multivariante:
El análisis de varianza (ANOVA) que se realizó a los datos se aplicaba para el modelo de parcelas divididas con DCA y no se empleó un MANOVA (análisis de varianza multivariante) debido a que en el artículo solo se evaluó una variable respuesta.
F.El método de comparaciones de medias después del Anova:
Según la lectura que se realizó al artículo , se puede ver que hace uso de un test de Tukey para evaluar los efectos individuales de cada factor, pero en la tabla de análisis de varianza se podía evidenciar que había interacción entre los factores, con un p.valor menor a 0.05. De forma que, el test de tukey se hizo para evaluar las diferencias significativas en las medias de cada factor de forma individual, y para evaluar las diferencias entre interacciones. Desde nuestro criterio no se debió aplicar el test de comparación de medias, para los factores por separado debido a que esa interacción significativa hace que los efectos individuales de un factor dependan de los niveles del otro factor.
G.La interacción de factores:
Según el p-valor observado en el ANOVA mostrado en los resultados del artículo existe una interacción entre el factor ubicado en la parcela principal y la subparcela, evidenciado con un p.valor menor al 0.05.
H.La presencia de bloques:
En este diseño, no hay presencia de bloques como tal, ya que se está usando el factor ambiente para bloquear al otro factor que es cultivar, característico de un diseño en parcelas divididas.
I.El balanceo o desbalanceo:
De los ensayos experimentales se obtuvo una base de datos con 90 observaciones, distribuidas en 6 ambientes. Esta base se construyó de tal forma que no existiesen datos faltantes, por lo tanto, es un diseño balanceado.
J.La definición clara de la unidad experimental:
La unidad experimental fue una parcela de 4 surcos de 6 m de largo, espaciados a 0,70 m entre sí.
K.Software utilizado y librería específica (en caso de ser R):
Las estimaciones de los modelos se llevaron a cabo mediante la utilización de los programas estadísticos InfoStat y Stata. Para el caso del análisis de varianza se utilizó InfoStat.