Crecimiento de pan.

En el libro Design and Analysis of Experiments with R de John Lawson se encuentra un ejemplo de un experimento que se realizó para determinar si el tiempo de reposo influía en el crecimiento de masas de panes.

A continuación se presentan los resultados de dicho experimento:

Lectura de la base de datos.

library(readxl)
pan <- read_excel("crecim_pan2.xlsx")

Convertir a tipo “factor” la variable explicativa.

pan$Tiempo <- factor(pan$Tiempo)
str(pan)
## Classes 'tbl_df', 'tbl' and 'data.frame':    12 obs. of  4 variables:
##  $ Idmasa: num  1 2 3 4 5 6 7 8 9 10 ...
##  $ Tiempo: Factor w/ 3 levels "35","40","45": 1 1 1 1 2 2 2 2 3 3 ...
##  $ Altura: num  4.5 5 5.5 6.75 6.5 6.5 10.5 9.5 9.75 8.75 ...
##  $ orden : num  8 12 11 1 9 10 7 5 4 6 ...

Gráfica exploratoria

library(ggplot2)
ggplot(pan, aes(Tiempo, Altura)) +
  geom_jitter(width = 0.1, height = 0) +
  labs(title = "Relación del tiempo de reposo con la altura alcanzada por el pan",
       x = "Tiempo (minutos)",
       y = "Altura (cm)") +
       scale_x_discrete(breaks = c("35","40","45")) +
       scale_y_continuous(breaks = seq(3,12)) +
       stat_summary(fun.y = "mean", geom = "point", col = "red", 
                    shape = 21, stroke = 2) +
       stat_summary(fun.y = "mean", geom = "line", col = "red",
                    linetype = 2, aes(group = 1))

Cálculo de estadísticos básicos

library(dplyr)
## 
## 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
pan %>%
  summarise(mediaA = mean(Altura), dsA = sd(Altura), numExp = n())
## # A tibble: 1 x 3
##   mediaA   dsA numExp
##    <dbl> <dbl>  <int>
## 1   7.33  1.97     12
pan %>% 
  group_by(Tiempo) %>% 
  summarise(mediaA = mean(Altura), dsA = sd(Altura),  numExp = n())
## # A tibble: 3 x 4
##   Tiempo mediaA   dsA numExp
##   <fct>   <dbl> <dbl>  <int>
## 1 35       5.44 0.966      4
## 2 40       8.25 2.06       4
## 3 45       8.31 1.36       4

Modelo.

Primera forma.

\[ y_{ij} = \mu_i + \epsilon_{ij} \] Donde: \[ \epsilon_{ij} \sim \mathcal{N}(0, \sigma^2) \quad i.i.d. \]

Segunda forma.

\[ \begin{eqnarray*} y_{ij} & = & \mu + \alpha_i + \epsilon_{ij} \\ \textrm{Dónde:} && \\ \epsilon_{ij} & \sim & \mathcal{N}(0, \sigma^2) \quad i.i.d. \end{eqnarray*} \]

Modelización.

modelo1aov <- aov(Altura ~ Tiempo, data = pan)
summary(modelo1aov)
##             Df Sum Sq Mean Sq F value Pr(>F)  
## Tiempo       2  21.57  10.786   4.602  0.042 *
## Residuals    9  21.09   2.344                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
modelo1lm <- lm(Altura ~ Tiempo, data = pan)
summary(modelo1lm)
## 
## Call:
## lm(formula = Altura ~ Tiempo, data = pan)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.812 -1.141  0.000  1.266  2.250 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   5.4375     0.7655   7.104 5.65e-05 ***
## Tiempo40      2.8125     1.0825   2.598   0.0288 *  
## Tiempo45      2.8750     1.0825   2.656   0.0262 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.531 on 9 degrees of freedom
## Multiple R-squared:  0.5056, Adjusted R-squared:  0.3958 
## F-statistic: 4.602 on 2 and 9 DF,  p-value: 0.042

Resultados de estimación de parámetros.

Promedios

model.tables(modelo1aov, type = "mean", se = TRUE)
## Tables of means
## Grand mean
##          
## 7.333333 
## 
##  Tiempo 
## Tiempo
##    35    40    45 
## 5.438 8.250 8.313 
## 
## Standard errors for differences of means
##         Tiempo
##          1.083
## replic.      4

Efectos

model.tables(modelo1aov, type = "effects", se = TRUE)
## Tables of effects
## 
##  Tiempo 
## Tiempo
##      35      40      45 
## -1.8958  0.9167  0.9792 
## 
## Standard errors of effects
##         Tiempo
##         0.7655
## replic.      4

Disgnósticos del modelo.

Diagnósticos gráficos preliminares.

plot(modelo1lm)

Normalidad de los residuales.

library(car)
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
id <- qqPlot(modelo1lm, id.n = 12, reps = 10000)

shapiro.test(rstudent(modelo1lm))
## 
##  Shapiro-Wilk normality test
## 
## data:  rstudent(modelo1lm)
## W = 0.94141, p-value = 0.5166

Prueba de igualdad de varianza

Si los residuales son normales.

bartlett.test(Altura ~ Tiempo, data = pan)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  Altura by Tiempo
## Bartlett's K-squared = 1.4732, df = 2, p-value = 0.4787

Si los residuales no son normales.

library(car)
leveneTest(modelo1lm)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  2.7992 0.1134
##        9
fligner.test(Altura ~ Tiempo, data = pan)
## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  Altura by Tiempo
## Fligner-Killeen:med chi-squared = 3.9103, df = 2, p-value = 0.1415
library(HH)
## Loading required package: lattice
## Loading required package: grid
## Loading required package: latticeExtra
## Loading required package: RColorBrewer
## 
## Attaching package: 'latticeExtra'
## The following object is masked from 'package:ggplot2':
## 
##     layer
## Loading required package: multcomp
## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
## 
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
## 
##     geyser
## Loading required package: gridExtra
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
## 
## Attaching package: 'HH'
## The following objects are masked from 'package:car':
## 
##     logit, vif
hov(Altura ~ Tiempo, data = pan)
## 
##  hov: Brown-Forsyth
## 
## data:  Altura
## F = 2.7992, df:Tiempo = 2, df:Residuals = 9, p-value = 0.1134
## alternative hypothesis: variances are not identical

Prueba de independencia

with(pan, plot(rstudent(modelo1lm) ~ orden, pch = 19,
               main = "Gráfica del orden de corrida",
               xlab = "Orden de corrida",
               ylab = "Residuales estandarizados",
               xlim = c(-1, 12),
               ylim = c(-2.1, 2.1)))
with(pan, text(orden, rstudent(modelo1lm), 
               labels  = paste(Idmasa, "t=", Tiempo), adj = 1,
               pos = 3, cex = 0.7))
abline(h = 0, col = "red", lty = 2)

library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
dwtest(modelo1lm, order.by = ~orden, data = pan, alternative = "two.sided")
## 
##  Durbin-Watson test
## 
## data:  modelo1lm
## DW = 1.5841, p-value = 0.3484
## alternative hypothesis: true autocorrelation is not 0

Comparaciones pareadas

Las compraraciones pareadas se dividen en dos clases:

Los métodos de Holm, Hochberg, Hommel y Bonferroni pertenecen a FWER. Intentan limitar la probabilidad de un falso positivo (cometer el error tipo I, rechazar Ho cuando en realidad es correcta), por lo tanto son conservadoras o muy estrictas.

Los métodos que pertenecen a FDR controlan el valor esperado de la proporción de falsos descubrimientos.

La mayoría de los métodos solo utilizan las probabilidades, sin embargo las pruebas de Tukey y Dunnet también utilizan la variabilidad de los datos. Las pruebas de Tukey y Dunnet se consideran métodos FWER.

No existe un consejo definitivo para determinar el ajuste a utilizar. En general se deberá escojer un método que le sea familiar a su audiencia de su campo de estudio. También se puede tener un poco de lógica a la hora de escoger algún tipo prueba, por ejemplo si es un estudio inicial, novedoso y preliminar se preferirían métodos que muestren significancia potencial para futuros estudios, si por el otro lado se tratan de estudios que exigen un algo grado de evidencia y que involucren vidas humanas como en caso médico se prefieren pruebas muy extrictas y conservadoras.

La prueba más liberal o más laxa es la que no hace ningún tipo de corrección y puede rechazar la hipótesis nula donde en realidad es cierta.

Por el otro lado la prueba más conservadora es la corrección de Bomferroni.

Entre estos dos extremos se encuentran diversas pruebas de comparación de medias.

La prueba de Tukey se considera dentro de las pruebas conservadoras.

Sin ajuste.

pairwise.t.test(x = pan$Altura, g = pan$Tiempo, p.adjust.method =  "none")
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  pan$Altura and pan$Tiempo 
## 
##    35    40   
## 40 0.029 -    
## 45 0.026 0.955
## 
## P value adjustment method: none

Con ajuste.

pairwise.t.test(x = pan$Altura, g = pan$Tiempo, p.adjust.method =  "bonferroni")
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  pan$Altura and pan$Tiempo 
## 
##    35    40   
## 40 0.086 -    
## 45 0.079 1.000
## 
## P value adjustment method: bonferroni
pairwise.t.test(x = pan$Altura, g = pan$Tiempo, p.adjust.method =  "hommel")
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  pan$Altura and pan$Tiempo 
## 
##    35    40   
## 40 0.058 -    
## 45 0.052 0.955
## 
## P value adjustment method: hommel
TukeyHSD(modelo1aov)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Altura ~ Tiempo, data = pan)
## 
## $Tiempo
##         diff        lwr      upr     p adj
## 40-35 2.8125 -0.2099347 5.834935 0.0676257
## 45-35 2.8750 -0.1474347 5.897435 0.0618204
## 45-40 0.0625 -2.9599347 3.084935 0.9981643