Introducción

Considere la siguiente pregunta: ¿Cúal de los nuevos fertilizantes presenta un mejor rendimiento considerando las condiciones de suelo?

De acuerdo a esto, lo que se busca es disminuir la variabilidad provocada por las variantes del tipo de suelo, buscando siempre el fertilizante que provoque un mejor resultado.

Entonces un diseño en Bloques al azar (DBA o RCBD) se usa cuando los experimentos tienen una fuente de variación conocida (Bloque) que podría afectar la respuesta, en general las fuentes en experimentos en terreno, son gradiente de sol, orientación, drenaje,etc.

Entonces cuando tienes varios tratamientos que quieres comparar, pero sabes que existe una condición extra (Bloque) que puede influir en el resultado, pero no es lo que quieres estudiar.

Entonces organizas el experimento en Bloques homogéneos, y dentro de cada Bloque aplicas todos los tratamientos al azar.

Así:

El Bloque controla la variabilidad “molesta”. Entonces el azar dentro de cada Bloque asegura imparcialidad.

Cuándo aplicarlo:

Cuando hay una fuente importante de variabilidad conocida que no puedes eliminar, pero puedes agrupar en Bloques (ejemplo: parcelas de Bloque con distinta fertilidad). Tambien cuando quieres aumentar la precisión de la comparación entre tratamientos.

Generación del diseño

El diseño en Bloques considera que existen al menos 2 factores o variables independientes, siendo una de estas más relevante que la otra.

el modelo a considera es

\(Y_{ij}= \mu + \tau_i + \beta_j + \epsilon\)

library(tidyverse)
library(emmeans)
library(rcompanion)
library(GGally)
library(nortest)
library(Analitica)
library(readxl)
library(writexl)

Contexto y objetivo

Objetivo: Evaluar el efecto de 5 tratamiento para en el rendimiento de un cultivo, Como contexto, se bloqueo en 3 condiciones, dado condiciones diferenciadas de suelo en eso 3 lugares, como mas cercanos a fuente de agua a más lejano.

Análisis de datos

Análisis Inicial

descripYG(datos,rendimiento,trat)
## Picking joint bandwidth of 96.6

##     Group  n      Mean Median         SD Kurtosis    Skewness         CV    Min
## 1 Control  9  123.9444  123.0   9.056919 1.841725  0.16284569 0.07307241  111.5
## 2   prod1  9 1607.5556 1522.0 257.494714 1.540550  0.28855539 0.16017780 1289.0
## 3   prod2  9 1308.8889 1299.0 179.045696 3.296203  0.82912577 0.13679213 1093.0
## 4   prod3  8  903.0000  898.0 189.573958 1.926825  0.32690600 0.20993794  659.0
## 5   prod4 10  897.5000  944.5 174.862009 1.804606 -0.40887100 0.19483232  618.0
## 6   prod5  9 1306.3333 1319.0 243.859591 1.699598  0.07856499 0.18667486 1009.0
##      Max     P25     P75    IQR
## 1  138.5  118.50  131.00  12.50
## 2 1969.0 1377.00 1897.00 520.00
## 3 1685.0 1212.00 1379.00 167.00
## 4 1204.0  748.00 1000.25 252.25
## 5 1131.0  759.75 1023.25 263.50
## 6 1688.0 1056.00 1490.00 434.00

Se observa una anomalía en los datos para el tratamiento Producto 2, al comparar con procedimiento anterior se ve el efecto ajuste grupo control que hacia disminuir el real rendimiento del cultivo en cada uno de los tratamientos.

analisis por replica (Bloque)

descripYG(datos,rendimiento,replica)
## Picking joint bandwidth of 221

##   Group  n     Mean Median       SD Kurtosis   Skewness        CV   Min  Max
## 1     1 18 1003.722 1068.5 505.3102 2.557086 -0.4761274 0.5034363 111.5 1897
## 2     2 18 1022.444 1130.0 513.7728 2.477988 -0.4710973 0.5024946 118.5 1902
## 3     3 18 1047.139 1058.5 543.7123 2.399304 -0.3615900 0.5192361 114.0 1969
##      P25     P75    IQR
## 1 780.75 1314.50 533.75
## 2 764.75 1357.50 592.75
## 3 777.75 1407.75 630.00

no se observa comportamiento anómalo desde la replcia, más allá de la sensibilidad provocada por el número de replicas o de bloques

Supuestos

#conversion de Variables Independientes en Factor
datos$trat<-as.factor(datos$trat)
datos$replica<-as.factor(datos$replica)

#verificacion
str(datos)
## tibble [54 × 3] (S3: tbl_df/tbl/data.frame)
##  $ trat       : Factor w/ 6 levels "Control","prod1",..: 3 3 3 3 3 3 3 3 3 1 ...
##  $ replica    : Factor w/ 3 levels "1","2","3": 1 1 1 2 2 2 3 3 3 1 ...
##  $ rendimiento: num [1:54] 1093 1301 1379 1212 1299 ...

Supuesto de normalidad

#modelo tentativo RBCD para la obtencion de residuales
modR<-lm(rendimiento~ replica+trat,data=datos)

#revisar aplicacion para ver si fue correctamente aplicado "df = al numero de grupos -1" en el caso de replicas 3 -1
anova(modR)
## Analysis of Variance Table
## 
## Response: rendimiento
##           Df   Sum Sq Mean Sq F value Pr(>F)    
## replica    2    17072    8536  0.2215 0.8022    
## trat       5 12081108 2416222 62.7024 <2e-16 ***
## Residuals 46  1772600   38535                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#obtener residuales para analisis
residuales<- data.frame(residual=modR$residuals)
descripYG(residuales,residual)

##    n          Mean    Median       SD Kurtosis  Skewness           CV       Min
## 1 54 -1.092387e-14 0.3716153 182.8805 2.305631 0.1867206 1.674136e+16 -320.4332
##        Max       P25      P75      IQR Fence_Low Fence_High
## 1 358.5668 -128.8962 91.43061 220.3268 -459.3864   421.9208
#Supuesto Normalidad

shapiro.test(residuales$residual)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuales$residual
## W = 0.96, p-value = 0.06887

La prueba entrega un p-value 0.06887 (mayor que alfa 0.05, algo algo cerca ), por tanto no se puede rechazar la hipótesis de que los residuales del modelo propuesto, presenten un comportamiento aproximadamente como Distribución de Probabilidad Normal.

Supuesto de Homocedásticidad

htrat<-Levene.Test(rendimiento~trat+replica, data = datos,center="mean")                  #
summary(htrat)
## 
## --- Homoscedasticity Test Summary ---
## 
## Method applied         : Levene (mean) - global by cells 
## F Statistic            : 1.517125 
## Degrees of freedom     : 17 (between), 36 (within)
## p-value                : 0.1437989 ns 
## Decision (alpha = 0.05): Homoscedastic (cells) 
## ----------------------------------------
## 
## >> ANOVA sobre desviaciones absolutas (descomposicion factorial)
## Metodo:                 ANOVA on |Y - cell_center| (mean), SS Type I (sequential)
## 
##       Term DF      F        p Signif   Interpretacion
##       trat  5 4.4823 0.002071     ** variance differs
##    replica  2 0.7618 0.472600     ns      no evidence
##  Residuals 46     NA       NA             no evidence
## ----------------------------------------

de acuerdo a la información se visualiza heterocedásticidad entre tratamientos. Lo que implica en anova una corrección de Welch, así mismo el procedimiento Post Hoc debe considerar la falta de homocedásticidad, de ser necesario, se recomienda Games Howell corregido para bloques.

Análisis del Anova respecto del rendimiento

Se considera un Diseño en Bloques al Azar, variable primaria al Tratamiento y Secundaria el Bloque (distancia con la fuente de agua).

# Permite modelar varianzas diferentes por tratamiento
modR<-lm(rendimiento~ replica+trat,data=datos)

#revisar aplicacion
anova(modR)
## Analysis of Variance Table
## 
## Response: rendimiento
##           Df   Sum Sq Mean Sq F value Pr(>F)    
## replica    2    17072    8536  0.2215 0.8022    
## trat       5 12081108 2416222 62.7024 <2e-16 ***
## Residuals 46  1772600   38535                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

\(H_0\) : los tratamiento NO provocan alguna diferencia en el rendimiento del cultivo

No se observa un efecto diferenciador provocado por la replica. Esto significa que No hay diferencias sistemáticas entre las replicas de los tratamiento aplicados al cultivo. Dado lo anterior se justifica la decisión de tomar Dosis como bloque, provocando una reducción de la variabilidad residual, mejorando la precisión del modelo.

En conclusión para el anova, No se puede aceptar la hipótesis \(H_0\), basado en que el p value es muy significativo (p < 0.000). Esto indica que habría una diferencia a nivel de tratamiento.
Entonces, existe al menos un tratamiento que provoca un mejor rendimiento independiente de la cercania con la fuente de agua.

Comparaciones entre tratamientos.

resGH<-pairwisePermutationTest(rendimiento ~ trat, 
                                        data = datos, 
                                        method = "fdr")
print(resGH)
##             Comparison    Stat   p.value  p.adjust
## 1  Control - prod1 = 0  -4.017 5.899e-05 0.0002636
## 2  Control - prod2 = 0  -4.042 5.307e-05 0.0002636
## 3  Control - prod3 = 0  -3.817  0.000135 0.0004050
## 4  Control - prod4 = 0   -4.05 5.119e-05 0.0002636
## 5  Control - prod5 = 0  -3.975 7.028e-05 0.0002636
## 6    prod1 - prod2 = 0   2.396   0.01656 0.0207000
## 7    prod1 - prod3 = 0   3.415 0.0006378 0.0013670
## 8    prod1 - prod4 = 0   3.669 0.0002437 0.0006092
## 9    prod1 - prod5 = 0   2.215   0.02674 0.0308500
## 10   prod2 - prod3 = 0   3.043  0.002344 0.0035550
## 11   prod2 - prod4 = 0    3.29  0.001003 0.0018810
## 12   prod2 - prod5 = 0 0.02612    0.9792 0.9792000
## 13   prod3 - prod4 = 0 0.06586    0.9475 0.9792000
## 14   prod3 - prod5 = 0   -2.79  0.005272 0.0071890
## 15   prod4 - prod5 = 0  -3.039   0.00237 0.0035550

de acuerdo a la test de Games Howell se observan 3 grupos de tratamientos, donde Producto 3 y Producto 4 presentan un comportamiento similar, así como Producto 2 y Producto 5

luego las diferencias de rendimiento del cultivo, se pueden observar en la siguiente gráfica.

resumen<-datos %>%
  group_by(trat) %>%
  summarise(Prom=mean(rendimiento),
            se = sd(rendimiento)/sqrt(n())
            )

letters_result <- cldList(p.adjust ~ Comparison,
                          data = resGH,
                          threshold = 0.05)

temp <- left_join(resumen, letters_result, 
                  by = c("trat" = "Group"))

temp <- temp %>% arrange(Prom)

temp$grupo_id <- as.numeric(factor(temp$Letter, levels = unique(temp$Letter)))

# Asignar nuevas letras en orden alfabético
temp$letra <- letters[temp$grupo_id]

# 6. Resultado final
resumen <- temp %>%
#  select(trat, Prom, letra) %>%
  arrange(Prom)
resumen$trat<-as.factor(resumen$trat)

ggplot(transform(resumen, trat = reorder(trat,
                                              Prom, 
                                              FUN = mean, na.rm = TRUE)),
       aes(x = trat, y = Prom, fill=trat)) +
  geom_bar(stat = "identity", width = 0.5) +
  geom_errorbar(aes(ymin = Prom - se, ymax = Prom + se), 
                width = 0.2, 
                color = "gray30") +
  ylab("rendimiento Promedio") +
  scale_y_continuous(limits = c(0, 2000), 
                     expand = expansion(mult = c(0, 0.05))) +
  ggtitle("rendimiento promedio en la col por tratamiento")+
   geom_text(aes(label = resumen$letra), 
            vjust = -2.5,
            size = 5) +
    theme_minimal() +
  theme(aspect.ratio = .8/1,  # Relación alto:ancho = 1:1.5
        plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 11),
        legend.position = "none")

Luego, el Producto 1, tendria un mayor efecto provocando un mayor rendimiento en el cultivo.