Retención de vitamina C en jugo de naranja congelado restituido

Se realizó un experimento en donde se midió el contenido de ácido ascórbico (vitamina C) en miligramos por litro de tres marcas de jugo de naranja en tres tiempos diferentes.

Las marcas fueron Rickfood, Sealed-Sweet y Minute Maid y se evaluó en los días cero (0), tres (3) y siete (7).

El tiempo se refiere al número de días a partir del momento que se exprimió el producto hasta que se descongeló para usarlo.

Se quiere saber si hay diferencia entre las marcas, el efecto del tiempo de congelación y la pérdida de vitamina C y si existe alguna interacción entre la marca y los días de congelación.

Se realizaron cuatro (4) réplicas para un total de 36 experimentos.

Cada jornada en donde se evaluó la cantidad de vitamina C, se tomó al azar el orden de medición de las 12 unidades del respectivo día.

Biblotecas de funciones utilizadas

library(ggplot2) ## Bibloteca para gráficas.
library(lmtest)  ## Biblioteca para algunas pruebas diagnósticas.
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(dplyr)   ## Biblioteca para manejo de base de datos.
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

Lectura de la base de datos

naranja<-read.csv2("naranja.csv")
naranja$dias <- factor(naranja$dias)

Análisis descriptivo

Cálculo de estadísticos básicos

sd(naranja$vitamC)
## [1] 3.895341
mean(naranja$vitamC)
## [1] 47.89444
summary(naranja$vitamC)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   40.40   44.00   48.30   47.89   49.65   56.00

Cálculo del coeficiente de variación

cv <- sd(naranja$vitamC)/mean(naranja$vitamC)*100
cv
## [1] 8.133179

Estadísticos por factores

## Cálculo de los promedio por un factor (dias)
medias1<-with(naranja, tapply(vitamC, list(dias), mean))
medias1
##        0        3        7 
## 51.25000 47.21667 45.21667
## Cálculo de los promedio por otro factor (tipoJugo)
medias2<-with(naranja, tapply(vitamC, list(tipoJugo), mean))
medias2
##  Minute Maid     Richfood Sealed-Sweet 
##     48.95000     48.10000     46.63333
## Cálculo de los promedios de cada combinación de los factores
medias<-with(naranja, tapply(vitamC, list(dias, tipoJugo), mean))
medias
##   Minute Maid Richfood Sealed-Sweet
## 0      52.475   50.775         50.5
## 3      48.200   48.650         44.8
## 7      46.175   44.875         44.6
## Cálculo de las desviaciones estándar de cada combinación de los factores
desvEstand<-with(naranja,tapply(vitamC, list(dias, tipoJugo), sd))
desvEstand
##   Minute Maid Richfood Sealed-Sweet
## 0   0.8057088 3.380705     3.729164
## 3   1.0708252 4.312385     2.771281
## 7   2.3157072 3.982775     3.174902
Otra forma de obtener los estadísticos por factor y combinado.
library(dplyr)

naranja %>% 
  group_by(dias) %>%
  summarise(promedio = mean(vitamC),
            desvEstand  = sd(vitamC))
## # A tibble: 3 x 3
##   dias  promedio desvEstand
##   <fct>    <dbl>      <dbl>
## 1 0         51.2       2.81
## 2 3         47.2       3.27
## 3 7         45.2       3.01
library(dplyr)

naranja %>% 
  group_by(tipoJugo) %>%
  summarise(promedio = mean(vitamC),
            desvEstand  = sd(vitamC))
## # A tibble: 3 x 3
##   tipoJugo     promedio desvEstand
##   <fct>           <dbl>      <dbl>
## 1 Minute Maid      49.0       3.08
## 2 Richfood         48.1       4.36
## 3 Sealed-Sweet     46.6       4.10
library(dplyr)

naranja %>% 
  group_by(dias, tipoJugo) %>%
  summarise(promedio = mean(vitamC),
            desvEstand  = sd(vitamC))
## # A tibble: 9 x 4
## # Groups:   dias [?]
##   dias  tipoJugo     promedio desvEstand
##   <fct> <fct>           <dbl>      <dbl>
## 1 0     Minute Maid      52.5      0.806
## 2 0     Richfood         50.8      3.38 
## 3 0     Sealed-Sweet     50.5      3.73 
## 4 3     Minute Maid      48.2      1.07 
## 5 3     Richfood         48.6      4.31 
## 6 3     Sealed-Sweet     44.8      2.77 
## 7 7     Minute Maid      46.2      2.32 
## 8 7     Richfood         44.9      3.98 
## 9 7     Sealed-Sweet     44.6      3.17

Análisis exploratorio de datos

## Gráfica utilizando el paquete de la gramática de gráficas.
library(ggplot2)
g1 <- ggplot(naranja, aes(dias, vitamC, col= tipoJugo))
g1 + geom_point(size = 3) +
     stat_summary(fun.y=mean, geom="line", size=1.2, aes(group = tipoJugo)) +
     ylab("Vitamina C (mg/l)") +
     xlab("Días") +
     ggtitle("Nivel de vitamina C vs días")

Modelo

El modelo para este caso es:

\[ y_{ijk} = \mu + \alpha_i + \beta_j + (\alpha \beta)_{ij} + \epsilon_{ijk} \\ i = 1, 2, 3. \quad j = 1, 2, 3. \quad \textrm{ y } \quad k = 1, 2, 3, 4. \\ \epsilon_{ijk} \sim \mathcal{N}(0, \sigma^2) i.i.d. \]

Hipótesis a probar

Factor marca de jugo

\[ H_0: \alpha_1 = \alpha_2 = \alpha_3 = 0 \\ H_1: \textrm{Alguna } \alpha \textrm{ diferente de 0} \]

Factor día

\[ H_0: \beta_1 = \beta_2 = \beta_3 = 0 \\ H_1: \textrm{Algún } \beta \textrm{ diferente de 0} \]

Interacción marca-días.

\[ H_0: (\alpha\beta)_{11} = \cdots = (\alpha\beta)_{33} = 0 \\ H_1: \textrm{Algún } (\alpha\beta)_{ij} \textrm{ diferente de 0} \]

Análisis de varianza

Modelo inicial

## Tabla anova para los dos factores y su interacción
m1A1<-aov(vitamC ~ dias + tipoJugo + dias:tipoJugo, data = naranja)
anova(m1A1)
## Analysis of Variance Table
## 
## Response: vitamC
##               Df  Sum Sq Mean Sq F value    Pr(>F)    
## dias           2 226.676 113.338 12.0411 0.0001827 ***
## tipoJugo       2  32.962  16.481  1.7510 0.1927491    
## dias:tipoJugo  4  17.301   4.325  0.4595 0.7646909    
## Residuals     27 254.140   9.413                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m1A2<-aov(vitamC ~ dias + tipoJugo, data = naranja)
anova(m1A2)
## Analysis of Variance Table
## 
## Response: vitamC
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## dias       2 226.676 113.338 12.9438 8.191e-05 ***
## tipoJugo   2  32.962  16.481  1.8822    0.1692    
## Residuals 31 271.441   8.756                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m1A3<-aov(vitamC ~ dias, data = naranja)
anova(m1A3)
## Analysis of Variance Table
## 
## Response: vitamC
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## dias       2 226.68 113.338  12.287 0.0001028 ***
## Residuals 33 304.40   9.224                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Promedios individuales por niveles del factor
model.tables(m1A3, type="mean", se = TRUE)
## Tables of means
## Grand mean
##          
## 47.89444 
## 
##  dias 
## dias
##     0     3     7 
## 51.25 47.22 45.22 
## 
## Standard errors for differences of means
##         dias
##         1.24
## replic.   12
## Cálculo de los efectos de cada nivel del factor
model.tables(m1A3, type="effects", se = TRUE)
## Tables of effects
## 
##  dias 
## dias
##      0      3      7 
##  3.356 -0.678 -2.678 
## 
## Standard errors of effects
##           dias
##         0.8768
## replic.     12
anovam2A<-anova(m1A3)
## Obtener sólo la desviación estándar de los errores
(dserror <- sqrt(anovam2A$"Mean Sq"[2]))
## [1] 3.03716

La estimación de la desviación estándar de los errores experimentales (\(\sigma\)) es de \(\hat{\sigma}=\) 3.04.

Diagnóstico del modelo

Prueba de normalidad

## Prueba de normalidad de los residuales
shapiro.test(residuals(m1A3))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(m1A3)
## W = 0.9597, p-value = 0.2104
## Prueba de normalidad de los residuales estandarizados
shapiro.test(rstandard(m1A3))
## 
##  Shapiro-Wilk normality test
## 
## data:  rstandard(m1A3)
## W = 0.9597, p-value = 0.2104
## Prueba de normalidad de los residuales estudentizados
shapiro.test(rstudent(m1A3))
## 
##  Shapiro-Wilk normality test
## 
## data:  rstudent(m1A3)
## W = 0.96397, p-value = 0.284

Prueba de igualdad de varianza

## Prueba de igualdad de varianza
## u homogeneidad de varianza o prueba de homocedasticidad.

## Prueba que exige normalidad, es la mejor bajo normalidad
bartlett.test(naranja$vitamC, naranja$dias)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  naranja$vitamC and naranja$dias
## Bartlett's K-squared = 0.24152, df = 2, p-value = 0.8862
library(lmtest)
## Prueba más robusta a desviaciones de la normal
bptest(naranja$vitamC ~ naranja$dias)
## 
##  studentized Breusch-Pagan test
## 
## data:  naranja$vitamC ~ naranja$dias
## BP = 0.55713, df = 2, p-value = 0.7569

Pruebas gráficas de aditividad y linealidad

## Pruebas gráficas
plot(m1A3)

Identificador de posible valores atípicos

residualesEstudentizados <- rstudent(m1A3)

idx <- which(abs(residualesEstudentizados) > qnorm(0.975))

plot(predict(m1A3), residualesEstudentizados, 
     ylim = c(-max(residualesEstudentizados), max(residualesEstudentizados)),
     pch = 19)
abline(h = c(qnorm(0.025), qnorm(0.975)), col = "red")
points(predict(m1A3)[idx], residualesEstudentizados[idx], pch = 19,
      col = "red")
text(predict(m1A3)[idx], residualesEstudentizados[idx], 
     labels = rownames(naranja)[idx], pos = 2)

Prueba de independencia

library(ggplot2)
g1 <- ggplot(naranja, aes(orden, residuals(m1A3)))
g1 + geom_point() +
     geom_smooth() +
     scale_x_continuous(name="Orden de ejecución", breaks = 1:12)
## `geom_smooth()` using method = 'loess'

library(ggplot2)
g1 <- ggplot(naranja, aes(orden, residuals(m1A3), col = tipoJugo))
g1 + geom_point() +
     geom_smooth() +
     scale_x_continuous(name="Orden de ejecución", breaks = 1:12)
## `geom_smooth()` using method = 'loess'

library(lmtest)
dwtest(vitamC ~ dias, order.by = ~ orden, data = naranja)
## 
##  Durbin-Watson test
## 
## data:  vitamC ~ dias
## DW = 1.8032, p-value = 0.3432
## alternative hypothesis: true autocorrelation is greater than 0

Resultados del modelo

Comparación de medias

compMed <- TukeyHSD(m1A3)
compMed
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = vitamC ~ dias, data = naranja)
## 
## $dias
##          diff       lwr        upr     p adj
## 3-0 -4.033333 -7.075831 -0.9908354 0.0072426
## 7-0 -6.033333 -9.075831 -2.9908354 0.0000796
## 7-3 -2.000000 -5.042498  1.0424980 0.2544508
plot(compMed)

Pronóstico del modelo

naranja$vitamCp <- predict(m1A3)

conf1 <- data.frame(predict(m1A3, interval="confidence"))
names(conf1) <- c("predM","infM","supM")
pred1 <- data.frame(predict(m1A3, interval="prediction"))
## Warning in predict.lm(m1A3, interval = "prediction"): predictions on current data refer to _future_ responses
names(pred1) <- c("predP","infP","supP")

naranja <- data.frame(naranja, conf1, pred1)
# Resumen de la predicción
aggregate(subset(naranja, select = c(vitamC:supP)),
          list(Días = naranja$dias), mean, na.rm=TRUE)
##   Días   vitamC  vitamCp    predM     infM     supM    predP     infP
## 1    0 51.25000 51.25000 51.25000 49.46623 53.03377 51.25000 44.81854
## 2    3 47.21667 47.21667 47.21667 45.43290 49.00043 47.21667 40.78520
## 3    7 45.21667 45.21667 45.21667 43.43290 47.00043 45.21667 38.78520
##       supP
## 1 57.68146
## 2 53.64813
## 3 51.64813
library(dplyr)
naranja %>% 
  group_by(dias) %>%
  summarise(promVitamC  = mean(vitamC),
            promVitamCp = mean(vitamCp),
            predM = mean(predM),
            infM  = mean(infM),
            supM  = mean(supM),
            predP = mean(predP),
            infP  = mean(infP),
            supP  = mean(supP))
## # A tibble: 3 x 9
##   dias  promVitamC promVitamCp predM  infM  supM predP  infP  supP
##   <fct>      <dbl>       <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0           51.2        51.2  51.2  49.5  53.0  51.2  44.8  57.7
## 2 3           47.2        47.2  47.2  45.4  49.0  47.2  40.8  53.6
## 3 7           45.2        45.2  45.2  43.4  47.0  45.2  38.8  51.6
library(ggplot2)
g1 <- ggplot(naranja, aes(dias, vitamC, col= tipoJugo))
g1 + geom_point(size = 3) +
     geom_line(aes(dias, vitamCp), col= "black", group = 1) +
     scale_y_continuous(breaks=seq(30,60,2)) +
     ylab("Vitamina C (mg/l)") +
     xlab("Días") +
     ggtitle("Nivel de vitamina C vs días.\nPredicción del promedio.")

library(ggplot2)
g1 <- ggplot(naranja, aes(dias, vitamC, col= tipoJugo))
limitesC <- aes(ymin = infM, ymax = supM)
g1 + geom_point(size = 3) +
     geom_line(aes(dias, vitamCp), col= "black", group = 1) +
     geom_errorbar(limitesC, width = 0.5, col = "black") +
     scale_y_continuous(breaks=seq(30,60,2)) +
     ylab("Vitamina C (mg/l)") +
     xlab("Días") +
     ggtitle("Nivel de vitamina C vs días.\nIntervalo de confianza para el promedio.")

library(ggplot2)
g1 <- ggplot(naranja, aes(dias, vitamC, col= tipoJugo))
limitesP <- aes(ymin = infP, ymax = supP)
g1 + geom_point(size = 3) +
     geom_line(aes(dias, vitamCp), col= "black", group = 1) +
     geom_errorbar(limitesP, width = 0.5, col = "black") +
     scale_y_continuous(breaks=seq(30,60,2))   +
     ylab("Vitamina C (mg/l)") +
     xlab("Días") +
     ggtitle("Nivel de vitamina C vs días.\nIntervalo de predicción para los datos.")