Ejercicio de rendimiento de almidon de yuca

Una empresa de alimentos está estudiando el efecto de 4 temperaturas sobre el rendimiento del almidón de yuca (en gramos). Cada evaluación tarda aproximadamente una hora y media, por lo que solamente se pueden realizar 4 evaluaciones por día. La idea general consiste en realizar una repetición completa del experimento en cada uno de los niveles del factor de control (factor bloque) y se realiza la aleatorización de los tratamientos dentro de cada bloque.

Temperatura Dia 1 Dia 2 Dia 3 Dia 4
30 14 14.1 14.5 14
35 13.9 13.8 14.2 14
40 14.1 14.2 14.4 13.9
45 13 13.1 13.5 13.2

Carga de datos

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.1     ✔ stringr   1.5.2
## ✔ ggplot2   4.0.0     ✔ tibble    3.3.0
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.1.0     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Yuca_desordenada <- read.csv("Yuca_desordenada.csv")
Yuca_desordenada
##   Temperatura Dia.1 Dia.2 Dia.3 Dia.4
## 1          30  14.0  14.1  14.5  14.0
## 2          35  13.9  13.8  14.2  14.0
## 3          40  14.1  14.2  14.4  13.9
## 4          45  13.0  13.1  13.5  13.2

Cambiando el punto . de los dias por un raya al piso _

colnames(Yuca_desordenada) <- c("Temperatura", "Dia_1", "Dia_2", "Dia_3", "Dia_4")
colnames(Yuca_desordenada)
## [1] "Temperatura" "Dia_1"       "Dia_2"       "Dia_3"       "Dia_4"

Estructurando los datos

yuca <- Yuca_desordenada %>% pivot_longer(cols = Dia_1:Dia_4, names_to = "Dia",
                                  values_to = "Rendimiento")
yuca
## # A tibble: 16 × 3
##    Temperatura Dia   Rendimiento
##          <int> <chr>       <dbl>
##  1          30 Dia_1        14  
##  2          30 Dia_2        14.1
##  3          30 Dia_3        14.5
##  4          30 Dia_4        14  
##  5          35 Dia_1        13.9
##  6          35 Dia_2        13.8
##  7          35 Dia_3        14.2
##  8          35 Dia_4        14  
##  9          40 Dia_1        14.1
## 10          40 Dia_2        14.2
## 11          40 Dia_3        14.4
## 12          40 Dia_4        13.9
## 13          45 Dia_1        13  
## 14          45 Dia_2        13.1
## 15          45 Dia_3        13.5
## 16          45 Dia_4        13.2

Analisis exploratorio de datos

yuca$Temperatura <- as.factor(yuca$Temperatura)
yuca$Dia <- as.factor(yuca$Dia)

# Resumen de los datos
str(yuca)
## tibble [16 × 3] (S3: tbl_df/tbl/data.frame)
##  $ Temperatura: Factor w/ 4 levels "30","35","40",..: 1 1 1 1 2 2 2 2 3 3 ...
##  $ Dia        : Factor w/ 4 levels "Dia_1","Dia_2",..: 1 2 3 4 1 2 3 4 1 2 ...
##  $ Rendimiento: num [1:16] 14 14.1 14.5 14 13.9 13.8 14.2 14 14.1 14.2 ...
summary(yuca)
##  Temperatura    Dia     Rendimiento   
##  30:4        Dia_1:4   Min.   :13.00  
##  35:4        Dia_2:4   1st Qu.:13.72  
##  40:4        Dia_3:4   Median :14.00  
##  45:4        Dia_4:4   Mean   :13.87  
##                        3rd Qu.:14.12  
##                        Max.   :14.50

Analisis grafico - Temperatura

bxp_temp <- yuca %>% ggplot(aes(x = Temperatura, y = Rendimiento, fill = Temperatura)) +
  geom_violin(width=0.8) +
  geom_boxplot(width=0.3, fill = "gray") +
  labs(x = "Temperatura (°C)", y  = "Rendimiento (%)", title = "Boxplot - Temperatura", caption = "UCEVA")

ggsave(plot = bxp_temp, "bxp_temp.png", dpi = 300, units = "in", height = 6, width = 8)

bxp_temp

Analisis grafico - Dias

bxp_dias <- yuca %>% ggplot(aes(x = Dia, y = Rendimiento, fill = Dia)) +
  geom_boxplot() +
  labs(x = "Dias", y  = "Rendimiento (%)", title = "Boxplot - Dias")

ggsave(plot = bxp_temp, "bxp_dias.png", dpi = 300, units = "in", height = 6, width = 8)

bxp_dias

Analisis grafico - Dias:Temperatura

bxp_dias_temp <- yuca %>% ggplot(aes(x = Dia, y = Rendimiento, col = Dia)) +
  facet_wrap(~Temperatura) +
  geom_point() +
  geom_line(aes(group = "Dia")) +
  theme(axis.text.x = element_text(angle=45, vjust=1, hjust=1)) +
  labs(x = NULL, y  = "Rendimiento (%)", title = "Cinetica de rendimiento - Dias:Temperatura")

ggsave(plot = bxp_dias_temp, "bxp_dias_temp.png", dpi = 300, units = "in", height = 6, width = 8)

bxp_dias_temp

Modelacion bajo un DBCA - ANOVA

modelo1 <- aov(Rendimiento ~ Temperatura + Dia, data = yuca)
summary(modelo1)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Temperatura  3 2.4669  0.8223   73.55 1.19e-06 ***
## Dia          3 0.4269  0.1423   12.73  0.00137 ** 
## Residuals    9 0.1006  0.0112                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Prueba post - ANOVA (Duncan)

Vamos a hacer la prueba para temperatura

library(agricolae)
compara_temp <- duncan.test(modelo1, "Temperatura")
compara_temp$groups %>% rownames_to_column("Temperatura")
##   Temperatura Rendimiento groups
## 1          30      14.150      a
## 2          40      14.150      a
## 3          35      13.975      a
## 4          45      13.200      b

Vamos a hacer la prueba para los dias (BLOQUES)

compara_dias <- duncan.test(modelo1, "Dia")
compara_dias$groups %>% rownames_to_column("Dia")
##     Dia Rendimiento groups
## 1 Dia_3      14.150      a
## 2 Dia_2      13.800      b
## 3 Dia_4      13.775      b
## 4 Dia_1      13.750      b

Prueba post - ANOVA (TUKEY)

Vamos a hacer la prueba para temperatura

compara_temp <- HSD.test(modelo1, "Temperatura")
compara_temp$groups %>% rownames_to_column("Temperatura")
##   Temperatura Rendimiento groups
## 1          30      14.150      a
## 2          40      14.150      a
## 3          35      13.975      a
## 4          45      13.200      b

Vamos a hacer la prueba para los dias (BLOQUES)

compara_dias <- HSD.test(modelo1, "Dia")
compara_dias$groups %>% rownames_to_column("Dia")
##     Dia Rendimiento groups
## 1 Dia_3      14.150      a
## 2 Dia_2      13.800      b
## 3 Dia_4      13.775      b
## 4 Dia_1      13.750      b

Resultado grafico dias

bar.group(compara_dias$groups, col = "blue", ylim=c(0,20))

Resultado grafico Tratamientos

bar.group(compara_temp$groups, col = "red",ylim=c(0,20))