1 Análisis estadístico

1.1 Carga de librerías

#-----------------------------------------------------------
# Instalación de paquetes
#-----------------------------------------------------------

paquetes <- c(
  "readxl",
  "dplyr",
  "ggplot2",
  "agricolae",
  "knitr"
)

instalar <- paquetes[!(paquetes %in% installed.packages()[,"Package"])]

if(length(instalar)>0){
  install.packages(instalar)
}

lapply(paquetes, library, character.only=TRUE)

## [[1]]
## [1] "readxl"    "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [7] "methods"   "base"     
## 
## [[2]]
## [1] "dplyr"     "readxl"    "stats"     "graphics"  "grDevices" "utils"    
## [7] "datasets"  "methods"   "base"     
## 
## [[3]]
##  [1] "ggplot2"   "dplyr"     "readxl"    "stats"     "graphics"  "grDevices"
##  [7] "utils"     "datasets"  "methods"   "base"     
## 
## [[4]]
##  [1] "agricolae" "ggplot2"   "dplyr"     "readxl"    "stats"     "graphics" 
##  [7] "grDevices" "utils"     "datasets"  "methods"   "base"     
## 
## [[5]]
##  [1] "knitr"     "agricolae" "ggplot2"   "dplyr"     "readxl"    "stats"    
##  [7] "graphics"  "grDevices" "utils"     "datasets"  "methods"   "base"

1.2 Razón de bloqueo

1.2.1 Gradiente de humedad

Figura 1. Cultivo de arveja durante el experimento.

1.3 Bloques realizados

Figura 1. Cultivo de arveja durante el experimento.

1.4 Importación de la base de datos

#-----------------------------------------------------------
# Lectura del archivo de Excel
#-----------------------------------------------------------

datos <- read_excel("Experimento arveja(1).xlsx")

head(datos)

## # A tibble: 6 × 6
##   Bloque Humedad Variedad      Planta Tratamiento  Peso
##    <dbl> <chr>   <chr>          <dbl> <chr>       <dbl>
## 1      1 Media   Pisum sativum      4 C           31.3 
## 2      1 Media   Pisum sativum      5 C           31.9 
## 3      1 Media   Pisum sativum      6 C            6.34
## 4      1 Media   Pisum sativum     10 T1          28.6 
## 5      1 Media   Pisum sativum     11 T1          54.6 
## 6      1 Media   Pisum sativum     12 T1          34.2

str(datos)

## tibble [36 × 6] (S3: tbl_df/tbl/data.frame)
##  $ Bloque     : num [1:36] 1 1 1 1 1 1 1 1 1 2 ...
##  $ Humedad    : chr [1:36] "Media" "Media" "Media" "Media" ...
##  $ Variedad   : chr [1:36] "Pisum sativum" "Pisum sativum" "Pisum sativum" "Pisum sativum" ...
##  $ Planta     : num [1:36] 4 5 6 10 11 12 16 17 18 7 ...
##  $ Tratamiento: chr [1:36] "C" "C" "C" "T1" ...
##  $ Peso       : num [1:36] 31.33 31.93 6.34 28.63 54.57 ...

1.5 Preparación de la base

#-----------------------------------------------------------
# Conversión de variables a factores
#-----------------------------------------------------------

datos$Bloque <- as.factor(datos$Bloque)

datos$Tratamiento <- as.factor(datos$Tratamiento)

datos$Variedad <- as.factor(datos$Variedad)

datos$Planta <- as.factor(datos$Planta)

summary(datos)

##  Bloque   Humedad                   Variedad      Planta   Tratamiento
##  1:9    Length:36          Pisum sativum:36   1      : 1   C :12      
##  2:9    Class :character                      2      : 1   T1:12      
##  3:9    Mode  :character                      3      : 1   T2:12      
##  4:9                                          4      : 1              
##                                               5      : 1              
##                                               6      : 1              
##                                               (Other):30              
##       Peso      
##  Min.   : 3.08  
##  1st Qu.:13.03  
##  Median :26.84  
##  Mean   :26.53  
##  3rd Qu.:33.52  
##  Max.   :77.07  
##

2 Estadística descriptiva

2.1 Estadísticos generales

summary(datos$Peso)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    3.08   13.03   26.84   26.53   33.52   77.07

cat("Número de observaciones:",length(datos$Peso),"\n")

## Número de observaciones: 36

cat("Media:",mean(datos$Peso),"\n")

## Media: 26.53028

cat("Mediana:",median(datos$Peso),"\n")

## Mediana: 26.845

cat("Desviación estándar:",sd(datos$Peso),"\n")

## Desviación estándar: 17.90754

cat("Coeficiente de variación:",
100*sd(datos$Peso)/mean(datos$Peso),"%\n")

## Coeficiente de variación: 67.4985 %

cat("Mínimo:",min(datos$Peso),"\n")

## Mínimo: 3.08

cat("Máximo:",max(datos$Peso),"\n")

## Máximo: 77.07

2.2 Estadísticos descriptivos por tratamiento

descriptivos_trat <- datos %>%

group_by(Tratamiento) %>%

summarise(

n=n(),

Media=mean(Peso),

Mediana=median(Peso),

DE=sd(Peso),

CV=100*sd(Peso)/mean(Peso),

Min=min(Peso),

Max=max(Peso)

)

kable(descriptivos_trat,
digits=2,
caption="Estadísticos descriptivos por tratamiento")

Estadísticos descriptivos por tratamiento
Tratamiento	n	Media	Mediana	DE	CV	Min	Max
C	12	24.30	19.65	18.68	76.86	6.34	77.07
T1	12	28.28	28.43	16.46	58.20	3.49	54.57
T2	12	27.01	26.99	19.79	73.26	3.08	72.41

2.3 Estadísticos descriptivos por bloque

descriptivos_bloque <- datos %>%

group_by(Bloque) %>%

summarise(

n=n(),

Media=mean(Peso),

Mediana=median(Peso),

DE=sd(Peso),

CV=100*sd(Peso)/mean(Peso),

Min=min(Peso),

Max=max(Peso)

)

kable(descriptivos_bloque,
digits=2,
caption="Estadísticos descriptivos por bloque")

Estadísticos descriptivos por bloque
Bloque	n	Media	Mediana	DE	CV	Min	Max
1	9	26.34	29.13	15.19	57.67	6.34	54.57
2	9	28.86	18.80	21.64	74.98	9.78	72.41
3	9	13.01	11.77	7.95	61.12	3.08	27.46
4	9	37.91	33.29	16.83	44.40	21.15	77.07

2.4 Boxplot por tratamiento

ggplot(datos,
aes(Tratamiento,
Peso,
fill=Tratamiento))+

geom_boxplot(alpha=0.75)+

theme_bw()+

labs(

title="Peso por tratamiento",

x="Tratamiento",

y="Peso"

)+

theme(legend.position="none")

2.5 Boxplot por bloque

ggplot(datos,
aes(Bloque,
Peso,
fill=Bloque))+

geom_boxplot(alpha=0.75)+

theme_bw()+

labs(

title="Peso por bloque",

x="Bloque",

y="Peso"

)+

theme(legend.position="none")

2.6 Histograma del rendimiento

ggplot(datos,
aes(Peso))+

geom_histogram(

bins=10,

color="black",

fill="steelblue"

)+

theme_bw()+

labs(

title="Distribución del peso",

x="Peso",

y="Frecuencia"

)

# Efectos de los tratamientos

## 10. Media general del experimento

#-----------------------------------------------------------
# Media general del experimento
#-----------------------------------------------------------

media_general <- mean(datos$Peso)

cat("Media general del experimento:", round(media_general,2))

## Media general del experimento: 26.53

2.7 Medias por tratamiento

#-----------------------------------------------------------
# Medias por tratamiento
#-----------------------------------------------------------

medias_trat <- datos %>%
  group_by(Tratamiento) %>%
  summarise(
    Media = mean(Peso)
  )

kable(
  medias_trat,
  digits = 2,
  caption = "Media del rendimiento por tratamiento"
)

Media del rendimiento por tratamiento
Tratamiento	Media
C	24.30
T1	28.28
T2	27.01

2.8 Estimación de los efectos de los tratamientos

#-----------------------------------------------------------
# Efecto = Media del tratamiento - Media general
#-----------------------------------------------------------

efectos_trat <- medias_trat %>%
  mutate(
    Efecto = Media - media_general
  )

kable(
  efectos_trat,
  digits = 3,
  caption = "Efectos estimados de los tratamientos"
)

Efectos estimados de los tratamientos
Tratamiento	Media	Efecto
C	24.303	-2.227
T1	28.276	1.746
T2	27.012	0.481

3 Efectos de los bloques

3.1 Medias por bloque

#-----------------------------------------------------------
# Medias por bloque
#-----------------------------------------------------------

medias_bloque <- datos %>%
  group_by(Bloque) %>%
  summarise(
    Media = mean(Peso)
  )

kable(
  medias_bloque,
  digits = 2,
  caption = "Media del rendimiento por bloque"
)

Media del rendimiento por bloque
Bloque	Media
1	26.34
2	28.86
3	13.01
4	37.91

3.2 Estimación de los efectos de los bloques

#-----------------------------------------------------------
# Efecto = Media del bloque - Media general
#-----------------------------------------------------------

efectos_bloque <- medias_bloque %>%
  mutate(
    Efecto = Media - media_general
  )

kable(
  efectos_bloque,
  digits = 3,
  caption = "Efectos estimados de los bloques"
)

Efectos estimados de los bloques
Bloque	Media	Efecto
1	26.342	-0.188
2	28.861	2.331
3	13.012	-13.518
4	37.906	11.375

3.3 Gráfico de medias por bloque

ggplot(efectos_bloque,
       aes(x = Bloque,
           y = Media,
           fill = Bloque)) +

  geom_col(width = 0.7) +

  geom_hline(yintercept = media_general,
             linetype = 2,
             colour = "red") +

  theme_bw() +

  labs(
    title = "Media del rendimiento por bloque",
    x = "Bloque",
    y = "Peso promedio"
  ) +

  theme(
    legend.position = "none"
  )

3.4 Interpretación de los efectos

cat("Interpretación:\n\n")

## Interpretación:

cat("• Un efecto positivo indica que el tratamiento produjo un rendimiento superior a la media general.\n")

## • Un efecto positivo indica que el tratamiento produjo un rendimiento superior a la media general.

cat("• Un efecto negativo indica que el tratamiento produjo un rendimiento inferior a la media general.\n")

## • Un efecto negativo indica que el tratamiento produjo un rendimiento inferior a la media general.

cat("• Para los bloques, un efecto positivo indica que el bloque presentó un rendimiento promedio superior al promedio general del experimento.\n")

## • Para los bloques, un efecto positivo indica que el bloque presentó un rendimiento promedio superior al promedio general del experimento.

#-----------------------------------------------------------
# Modelo DBCA
#-----------------------------------------------------------

modelo <- aov(Peso ~ Tratamiento + Bloque,
              data = datos)

3.5 Tabla ANOVA

#-----------------------------------------------------------
# Tabla ANOVA
#-----------------------------------------------------------

anova_resultados <- summary(modelo)

anova_resultados

##             Df Sum Sq Mean Sq F value Pr(>F)  
## Tratamiento  2     99    49.4   0.179 0.8367  
## Bloque       3   2858   952.8   3.458 0.0286 *
## Residuals   30   8267   275.6                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.6 Interpretación del ANOVA

#-----------------------------------------------------------
# Interpretación automática
#-----------------------------------------------------------

p_trat <- anova_resultados[[1]]$`Pr(>F)`[1]

if(p_trat < 0.05){

cat("Conclusión:\n\n")

cat("Se rechaza H0.\n")

cat("Existen diferencias significativas entre tratamientos.\n")

cat("Se procede con la prueba de Tukey.\n")

}else{

cat("Conclusión:\n\n")

cat("No se rechaza H0.\n")

cat("No existen diferencias significativas entre tratamientos.\n")

}

## Conclusión:
## 
## No se rechaza H0.
## No existen diferencias significativas entre tratamientos.

4 Comparación múltiple de medias (Tukey)

4.1 Prueba de Tukey

#-----------------------------------------------------------
# Comparación múltiple
#-----------------------------------------------------------

tukey <- TukeyHSD(modelo,
                  which = "Tratamiento")

tukey

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Peso ~ Tratamiento + Bloque, data = datos)
## 
## $Tratamiento
##            diff       lwr      upr     p adj
## T1-C   3.972500 -12.73416 20.67916 0.8285165
## T2-C   2.708333 -13.99833 19.41500 0.9159729
## T2-T1 -1.264167 -17.97083 15.44250 0.9810110

4.2 Gráfico de Tukey

plot(tukey,
     las = 1,
     col = "steelblue")

4.3 Tabla de comparaciones

comparaciones <- as.data.frame(tukey$Tratamiento)

comparaciones

##            diff       lwr      upr     p adj
## T1-C   3.972500 -12.73416 20.67916 0.8285165
## T2-C   2.708333 -13.99833 19.41500 0.9159729
## T2-T1 -1.264167 -17.97083 15.44250 0.9810110

4.4 Medias de los tratamientos

medias_trat

## # A tibble: 3 × 2
##   Tratamiento Media
##   <fct>       <dbl>
## 1 C            24.3
## 2 T1           28.3
## 3 T2           27.0

4.5 Resumen de Tukey

cat("Interpretación:\n\n")

## Interpretación:

cat("Cada fila representa la comparación entre dos tratamientos.\n\n")

## Cada fila representa la comparación entre dos tratamientos.

cat("La columna 'diff' corresponde a la diferencia entre medias.\n\n")

## La columna 'diff' corresponde a la diferencia entre medias.

cat("La columna 'p adj' representa el valor-p ajustado mediante el método de Tukey.\n\n")

## La columna 'p adj' representa el valor-p ajustado mediante el método de Tukey.

cat("Cuando p adj < 0.05, las medias de esos tratamientos son significativamente diferentes.")

## Cuando p adj < 0.05, las medias de esos tratamientos son significativamente diferentes.


``` r
# Verificación de supuestos

## Extracción de residuos

#-----------------------------------------------------------
# Residuos y valores ajustados
#-----------------------------------------------------------

residuos <- residuals(modelo)

ajustados <- fitted(modelo)

4.6 Prueba de normalidad (Shapiro-Wilk)

Hipótesis

H₀: Los residuos siguen una distribución normal.
H₁: Los residuos no siguen una distribución normal.

#-----------------------------------------------------------
# Prueba de Shapiro-Wilk
#-----------------------------------------------------------

shapiro.test(residuos)

## 
##  Shapiro-Wilk normality test
## 
## data:  residuos
## W = 0.8957, p-value = 0.002608

4.7 Homogeneidad de varianzas (Bartlett)

Hipótesis

H₀: Las varianzas son homogéneas.
H₁: Al menos una varianza difiere.

#-----------------------------------------------------------
# Bartlett
#-----------------------------------------------------------

bartlett.test(Peso ~ Tratamiento,
              data = datos)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Peso by Tratamiento
## Bartlett's K-squared = 0.36666, df = 2, p-value = 0.8325

5 Conclusiones

5.1 Conclusión general del análisis

cat("=========================================\n")

## =========================================

cat(" RESUMEN DEL ANÁLISIS ESTADÍSTICO\n")

##  RESUMEN DEL ANÁLISIS ESTADÍSTICO

cat("=========================================\n\n")

## =========================================

#-----------------------------------------
# ANOVA
#-----------------------------------------

p_trat <- anova_resultados[[1]]$`Pr(>F)`[1]

if(p_trat < 0.05){

cat("ANOVA:\n")

cat("- Se detectaron diferencias significativas entre tratamientos.\n")

cat("- Se recomienda interpretar la prueba de Tukey.\n\n")

}else{

cat("ANOVA:\n")

cat("- No se detectaron diferencias significativas entre tratamientos.\n\n")

}

## ANOVA:
## - No se detectaron diferencias significativas entre tratamientos.

#-----------------------------------------
# Shapiro
#-----------------------------------------

p_shapiro <- shapiro.test(residuos)$p.value

if(p_shapiro > 0.05){

cat("Normalidad:\n")

cat("- Los residuos presentan distribución normal.\n\n")

}else{

cat("Normalidad:\n")

cat("- Los residuos no presentan distribución normal.\n\n")

}

## Normalidad:
## - Los residuos no presentan distribución normal.

#-----------------------------------------
# Bartlett
#-----------------------------------------

p_bartlett <- bartlett.test(Peso ~ Tratamiento,
                            data = datos)$p.value

if(p_bartlett > 0.05){

cat("Homogeneidad:\n")

cat("- Se cumple el supuesto de igualdad de varianzas.\n\n")

}else{

cat("Homogeneidad:\n")

cat("- No se cumple el supuesto de igualdad de varianzas.\n\n")

}

## Homogeneidad:
## - Se cumple el supuesto de igualdad de varianzas.

Efecto de diferentes dosis de caolín (Al₂Si₂O₅(OH)₄) sobre el rendimiento de grano por planta de arveja (Pisum Sativum)

Diseño en Bloques Completos al Azar (DBCA)

July Acosta, Santiago Diaz, José Rincón, Paula Valbuena

1 Análisis estadístico

1.1 Carga de librerías

1.2 Razón de bloqueo

1.2.1 Gradiente de humedad

1.3 Bloques realizados

1.4 Importación de la base de datos

1.5 Preparación de la base

2 Estadística descriptiva

2.1 Estadísticos generales

2.2 Estadísticos descriptivos por tratamiento

2.3 Estadísticos descriptivos por bloque

2.4 Boxplot por tratamiento

2.5 Boxplot por bloque

2.6 Histograma del rendimiento

2.7 Medias por tratamiento

2.8 Estimación de los efectos de los tratamientos

3 Efectos de los bloques

3.1 Medias por bloque

3.2 Estimación de los efectos de los bloques

3.3 Gráfico de medias por bloque

3.4 Interpretación de los efectos

3.5 Tabla ANOVA

3.6 Interpretación del ANOVA

4 Comparación múltiple de medias (Tukey)

4.1 Prueba de Tukey

4.2 Gráfico de Tukey

4.3 Tabla de comparaciones

4.4 Medias de los tratamientos

4.5 Resumen de Tukey

4.6 Prueba de normalidad (Shapiro-Wilk)

4.7 Homogeneidad de varianzas (Bartlett)

5 Conclusiones

5.1 Conclusión general del análisis