Diseño por bloques completamente aleatorizado (DBCA)

Paquetes necesario para la lectura de archivos excel

library(readxl)

## Warning: package 'readxl' was built under R version 4.4.2

a) Identificación del Factor de Tratamiento y Factor de Bloque

Factor de Tratamiento:
- Aceite: Este es el factor de tratamiento en el diseño experimental. Se refiere a los diferentes tipos de aceite utilizados en el experimento. Cada tipo de aceite representa un nivel del factor de tratamiento, y se espera que tenga un efecto en la variable de respuesta (por ejemplo, el ahorro de combustible).
Factor de Bloque:
- Camión: Este es el factor de bloque en el diseño experimental. Se refiere a los diferentes camiones utilizados en el experimento. El uso de camiones como factor de bloque ayuda a controlar la variabilidad que puede surgir de las diferencias entre los camiones, permitiendo así una comparación más precisa de los efectos del tratamiento (aceite) sobre la variable de respuesta.

Resumen

Factor de Tratamiento: Aceite
Factor de Bloque: Camión

library(openxlsx)
dat1<- read.xlsx("C:/rstudio/segundodbca.xlsx")
print(dat1)

##    Aceite tratamiento     y
## 1       1           1   0.5
## 2       1           2 0.634
## 3       1           3 0.487
## 4       1           4 0.329
## 5       1           5 0.512
## 6       2           1 0.535
## 7       2           2 0.675
## 8       2           3  0.52
## 9       2           4 0.435
## 10      2           5  0.54
## 11      3           1 0.513
## 12      3           2 0.595
## 13      3           3 0.488
## 14      3           4   0.4
## 15      3           5  0.51

# Verificar el tipo de datos de 'y'
str(dat1$y)

##  chr [1:15] "0.5" "0.634" "0.487" "0.329" "0.512" "0.535" "0.675" "0.52" ...

dat1$Aceite <- as.factor(dat1$Aceite)  # Si Aceite está presente y es categórica
dat1$tratamiento <- as.factor(dat1$tratamiento)

dat1$y <- as.numeric(dat1$y)
# Verificar si hay valores faltantes (NA)
na_y <- sum(is.na(dat1$y))
na_tratamiento <- sum(is.na(dat1$tratamiento))
cat("Número de valores NA en 'y':", na_y, "\n")

## Número de valores NA en 'y': 0

cat("Número de valores NA en 'tratamiento':", na_tratamiento, "\n")

## Número de valores NA en 'tratamiento': 0

if (na_y > 0 || na_tratamiento > 0) {
  print(dat1[!complete.cases(dat1), ])
  
  # Opcional: Eliminar filas con NA
  dat1 <- na.omit(dat1)
  cat("Se eliminaron las filas con NA.\n")
}

boxplot(y ~ tratamiento, data = dat1, 
        xlab = "Tratamiento",  
        ylab = "Ahorro de Combustible", 
        main = "Boxplot del Ahorro de Combustible por Tratamiento",
        las = 1, 
        cex.axis = 1.2, 
        cex.lab = 1.2)

Análisis del Box Plot

Comparación de Tratamientos

El box plot muestra la distribución del ahorro de combustible para cada tratamiento.
Se observan diferencias significativas entre los tratamientos en términos de la mediana y la dispersión.
Algunos tratamientos (como el 4) tienen una mediana notablemente más baja, lo que podría indicar que generan un menor ahorro de combustible.
El tratamiento 2 parece ser el más efectivo, con una mediana más alta y menor dispersión.

Dispersión y Outliers

La dispersión (tamaño de los cuadros y líneas) varía entre tratamientos.
No se observan outliers evidentes, lo que sugiere que no hay valores atípicos extremos.

DBCA

ANOVA

b) Hipótesis y modelo estadístico

El modelo para un DBCA es:

\[ Y_{ij} = \mu + \tau_i + \beta_j + \epsilon_{ij} \]

Donde:

$\tau_i$: efecto del i-ésimo tratamiento (Aceite).
$\beta_j$: efecto del j-ésimo bloque (Camión).
$\epsilon_{ij}$: error aleatorio.

Hipótesis

Hipótesis nula ($H_0$): No hay diferencia en el efecto de los diferentes tratamientos de aceite sobre la variable de respuesta.
Hipótesis alternativa ($H_a$): Al menos uno de los tratamientos de aceite tiene un efecto diferente sobre la variable de respuesta.

c) ANOVA

fit1 = lm(y ~ tratamiento  + Aceite, data=dat1)
doe1=anova(fit1)
doe1

Interpretación del ANOVA

La tabla ANOVA te permite determinar si las diferencias observadas entre los tratamientos y los bloques (Aceite) son estadísticamente significativas.

Tratamiento

Tiene un valor F alto (43.63) y un p-valor de 1.78e-05 (extremadamente bajo), indicando que existen diferencias significativas entre los tratamientos en términos del ahorro de combustible.

Aceite (Bloques)

El factor Aceite también es significativo (p-valor = 0.0223), lo que significa que las variaciones entre bloques también contribuyen significativamente al modelo.

Error Residual

La suma de cuadrados residuales es baja, lo que indica que el modelo ajusta bien los datos.

Conclusiones

Existe evidencia estadística fuerte para afirmar que tanto el tratamiento como el tipo de aceite tienen un efecto significativo sobre el ahorro de combustible.
Basándote en el box plot y ANOVA, se recomendaría analizar más de cerca los tratamientos específicos para identificar cuál es el más eficiente y realizar comparaciones múltiples para confirmar las diferencias.

Hipótesis

$H_0$: No hay diferencias significativas entre los tratamientos.
$H_1$: Al menos un tratamiento es diferente.

La tabla ANOVA muestra que el p-valor para el efecto del tratamiento es 1.781e-05, que es mucho menor que 0.05. Esto indica que hay evidencia estadística suficiente para rechazar la hipótesis nula $H_0$, que plantea que no hay diferencias significativas entre los tratamientos.

Por lo tanto, se concluye que existen diferencias significativas en el ahorro de combustible entre los diferentes tratamientos (aceites).

Los grados de libertad para cada uno de las sumas de los cuadrados esta contenido en anova(fit1)$Df

# grados de libertad
gl_anova=anova(fit1)$Df

Los grados de libertad gl para los tratamientos, para los bloques y el error son respectivamente

# tratamientos
gl.trat=gl_anova[1]
# bloques
gl.bl=gl_anova[2]
# error
gl.error=gl_anova[3]

En anova(fit1)$Sum se almacenan la suma de los cuadrados, SC.

sc_anova=anova(fit1)$Sum

Luego, la SC para los tratamientos, para los bloques y el error son respectivamente

sctr=sc_anova[1];sctr

## [1] 0.09209973

scbl=sc_anova[2];scbl

## [1] 0.006705733

sce=sc_anova[2];sce

## [1] 0.006705733

En anova(fit1)$Mean se almacenan los cuadrados medios, CM

cm_anova=anova(fit1)$Mean

Luego, los CM para los tratamientos, para los bloques y el error son respectivamente

cmtrat=cm_anova[1]
cmbl=cm_anova[2]
cme=cm_anova[3]

Verificación del ANOVA

doe1 = anova(fit1)
doe1

COMPARACIONES MULTIPLES

Metodo Fisher-LSD

library(agricolae)

## Warning: package 'agricolae' was built under R version 4.4.2

#LSD.test( fit1, "tratamiento", DFerror=gl.error, MSerror=cme, p.adj = "none")
lsd_dbc1=LSD.test(dat1$y, dat1$tratamiento, DFerror=gl.error, MSerror=cme, p.adj = "none")
lsd_dbc1

## $statistics
##        MSerror Df      Mean       CV  t.value        LSD
##   0.0005277833  8 0.5115333 4.491112 2.306004 0.04325559
## 
## $parameters
##         test p.ajusted           name.t ntr alpha
##   Fisher-LSD      none dat1$tratamiento   5  0.05
## 
## $means
##      dat1$y        std r         se       LCL       UCL   Min   Max    Q25
## 1 0.5160000 0.01769181 3 0.01326378 0.4854137 0.5465863 0.500 0.535 0.5065
## 2 0.6346667 0.04000417 3 0.01326378 0.6040803 0.6652530 0.595 0.675 0.6145
## 3 0.4983333 0.01877054 3 0.01326378 0.4677470 0.5289197 0.487 0.520 0.4875
## 4 0.3880000 0.05400926 3 0.01326378 0.3574137 0.4185863 0.329 0.435 0.3645
## 5 0.5206667 0.01677299 3 0.01326378 0.4900803 0.5512530 0.510 0.540 0.5110
##     Q50    Q75
## 1 0.513 0.5240
## 2 0.634 0.6545
## 3 0.488 0.5040
## 4 0.400 0.4175
## 5 0.512 0.5260
## 
## $comparison
## NULL
## 
## $groups
##      dat1$y groups
## 2 0.6346667      a
## 5 0.5206667      b
## 1 0.5160000      b
## 3 0.4983333      b
## 4 0.3880000      c
## 
## attr(,"class")
## [1] "group"

# grafico
plot( lsd_dbc1, las=1,cex=1.5, cex.axis=1.3)

Metodo de Tukey

TukeyHSD(aov(fit1), "tratamiento")

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = fit1)
## 
## $tratamiento
##             diff         lwr         upr     p adj
## 2-1  0.118666667  0.05386311  0.18347022 0.0015207
## 3-1 -0.017666667 -0.08247022  0.04713689 0.8729154
## 4-1 -0.128000000 -0.19280355 -0.06319645 0.0009128
## 5-1  0.004666667 -0.06013689  0.06947022 0.9989608
## 3-2 -0.136333333 -0.20113689 -0.07152978 0.0005914
## 4-2 -0.246666667 -0.31147022 -0.18186311 0.0000076
## 5-2 -0.114000000 -0.17880355 -0.04919645 0.0019827
## 4-3 -0.110333333 -0.17513689 -0.04552978 0.0024541
## 5-3  0.022333333 -0.04247022  0.08713689 0.7569059
## 5-4  0.132666667  0.06786311  0.19747022 0.0007141

Grafico para las comparaciones múltiples

plot(TukeyHSD(aov(fit1), "tratamiento"),las=1, col=1)

Resultados del Análisis ANOVA

Los resultados del análisis ANOVA confirman que existen diferencias significativas entre los aceites ($p = 0.022$), lo que indica que no todos los aceites tienen el mismo impacto en la variable de respuesta ($y$). Además, se encontraron diferencias altamente significativas entre los tratamientos ($p = 1.78 \times 10^{-5}$).

Las pruebas post hoc (LSD y Tukey) confirman que el Tratamiento 2 (Aceite 2) presenta un rendimiento significativamente superior en comparación con los demás aceites. Esto indica que la elección del tipo de aceite influye de manera importante en la eficiencia.

En contraste, el Tratamiento 4 (Aceite 4) presenta el rendimiento más bajo, mientras que los tratamientos correspondientes a Aceites 1, 3 y 5 no mostraron diferencias significativas entre sí.

VERIFICACION DE HIPOTESIS

Residuos

Los residuos son obtenidos

# residuos
resi_fit1=fit1$residuals

Los valores ajustados

# valores ajustados
yhat_fit1=fit1$fitted.values

Los valores observados, ajustados y residuos

# datos observados, ajuste residuos
cbind(dat1$y, yhat_fit1, resi_fit1)

##          yhat_fit1     resi_fit1
## 1  0.500 0.4968667  3.133333e-03
## 2  0.634 0.6155333  1.846667e-02
## 3  0.487 0.4792000  7.800000e-03
## 4  0.329 0.3688667 -3.986667e-02
## 5  0.512 0.5015333  1.046667e-02
## 6  0.535 0.5454667 -1.046667e-02
## 7  0.675 0.6641333  1.086667e-02
## 8  0.520 0.5278000 -7.800000e-03
## 9  0.435 0.4174667  1.753333e-02
## 10 0.540 0.5501333 -1.013333e-02
## 11 0.513 0.5056667  7.333333e-03
## 12 0.595 0.6243333 -2.933333e-02
## 13 0.488 0.4880000 -1.550409e-17
## 14 0.400 0.3776667  2.233333e-02
## 15 0.510 0.5103333 -3.333333e-04

Normalidad de los errores

El gráfico de normalidad para los errores. Se observa que los errores siguen aproximadamente una distribución normal

# grafica de probabilidad normal 
plot(fit1,which=2, col="red", pch=19 )

Se puede usar paquetes para verificar la normalidad de los errores tales como

paquete car
paquete performance

Usando el paquete car

car::qqPlot(fit1$residuals, distribution="norm", id=FALSE, ylab="residuales",
            las=1)

Usando el paquete performance

library(performance)
norm_fit1=check_normality(fit1)
plot(norm_fit1, type="qq")

Se puede usar pruebas de hipótesis para verificar la normalidad de los errores. En R se puede usar los test de Shapiro-Wilk(SW) y de Kolmogorov-Smirnov(KS).

Para la prueba de hipótesis se tiene
$H_0 :$ Los errores son normales
$H_1 :$ Los errores no son normales

e) Verificación de Supuestos del Modelo

La verificación de los supuestos del modelo es un paso crucial en el análisis de ANOVA, ya que garantiza la validez de los resultados obtenidos. Los principales supuestos a verificar son:

Normalidad de los Residuos:
- Los residuos del modelo deben seguir una distribución normal. Esto se puede evaluar mediante:
  - Prueba de Shapiro-Wilk: Para determinar si los residuos se distribuyen normalmente.
Homogeneidad de Varianzas:
- Las varianzas de los diferentes grupos deben ser aproximadamente iguales. Esto se puede comprobar mediante:
  - Prueba de Levene: Para evaluar si las varianzas son iguales entre los grupos.
Independencia de las Observaciones:
- Las observaciones deben ser independientes entre sí. Este supuesto se asegura generalmente mediante un diseño experimental adecuado.

Procedimientos de Verificación

Normalidad: Realizar la prueba de Shapiro-Wilk.
Homogeneidad de Varianzas: Aplicar la prueba de Levene.
Independencia: Asegurarse de que el diseño experimental no introduzca sesgos que afecten la independencia de las observaciones.

shapiro.test(fit1$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  fit1$residuals
## W = 0.91855, p-value = 0.183

Prueba de Normalidad de Shapiro-Wilk

Se realizó la prueba de normalidad de Shapiro-Wilk sobre los residuos del modelo para evaluar si estos siguen una distribución normal. Los resultados de la prueba son los siguientes:

Estadístico W: 0.91855
p-valor: 0.183

Interpretación

Dado que el p-valor (0.183) es mayor que el nivel de significancia comúnmente utilizado (0.05), no se rechaza la hipótesis nula de que los residuos siguen una distribución normal. Esto sugiere que los residuos del modelo son aproximadamente normales, lo que es un requisito para la validez del ANOVA.

ks.test(fit1$residuals,"pnorm")

## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  fit1$residuals
## D = 0.49109, p-value = 0.0007347
## alternative hypothesis: two-sided

Prueba de Kolmogorov-Smirnov

Se realizó la prueba de Kolmogorov-Smirnov para evaluar si los residuos del modelo siguen una distribución normal. Los resultados de la prueba son los siguientes:

Estadístico D: 0.49109
p-valor: 0.0007347
Hipótesis alternativa: dos lados

Interpretación

Dado que el p-valor (0.0007347) es significativamente menor que el nivel de significancia comúnmente utilizado (0.05), se rechaza la hipótesis nula de que los residuos siguen una distribución normal. Esto sugiere que los residuos del modelo no se distribuyen normalmente, lo que puede afectar la validez de los resultados del análisis ANOVA y otros métodos estadísticos que asumen normalidad.

Conclusiones

La falta de normalidad en los residuos puede indicar la necesidad de considerar transformaciones de los datos o el uso de métodos estadísticos que no requieran este supuesto. ## Evaluación de la Normalidad de los Residuos

Se realizaron dos pruebas de normalidad para evaluar si los residuos del modelo siguen una distribución normal:

Prueba de Kolmogorov-Smirnov: Detectó desviaciones significativas de la normalidad (p = 0.0007347).
Prueba de Shapiro-Wilk: No encontró evidencias significativas para rechazar la normalidad (p = 0.183).

Dado que el tamaño de la muestra es pequeño (n = 15), se da mayor peso a los resultados de la prueba de Shapiro-Wilk, que es más robusta en estas condiciones. Por lo tanto, se concluye que los residuos siguen una distribución normal. Sin embargo, esta conclusión puede complementarse con gráficos de diagnóstico para confirmar visualmente la normalidad.

Independencia de los errores

plot(fit1$residuals, col=2,cex=1.4, las=1, pch=19, xlab = "orden", ylab = "residuos")

Evaluación de la Independencia de los Residuos

La gráfica indica que los residuos no muestran evidencia de dependencia y parecen comportarse de manera independiente. Por lo tanto, podemos concluir que el supuesto de independencia de los residuos se cumple.

Homogeneidad de varianzas (Homocedasticidad)

#residuals vs yhats
plot(fit1,which=1)

Pruebas estadisticas para la homocedasticidad

Para hacer las pruebas estadisticas para verificar la homogeneidad de las varianzas se puede realizarlo usando las pruebas de

Bartlett: asumiendo normalidad
Levene: Es robusto y no presupone la normalidad de los datos.

En R se tiene que usar algunos paquetes para este fin tales como car, lawstat, EnvStats, ‘vartest enre otros

Prueba de Bartlett

# PRUEBA DE BARLETT
bartlett.test(y~tratamiento, data=dat1)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  y by tratamiento
## Bartlett's K-squared = 3.9598, df = 4, p-value = 0.4115

Resultados de la Prueba de Bartlett

Estadístico de Bartlett (K²): 3.9598
Grados de libertad (df): 4
p-valor: 0.4115

Dado que el p-valor es mayor que el nivel de significancia común (α = 0.05), no se rechaza la hipótesis nula de que las varianzas son iguales.

usando el paquete EnvStats

EnvStats::varGroupTest( y~tratamiento, data = dat1, test = "Bartlett")

## 
## Results of Hypothesis Test
## --------------------------
## 
## Null Hypothesis:                 Ratio of each pair of variances = 1
## 
## Alternative Hypothesis:          At least one variance differs
## 
## Test Name:                       Bartlett's Test for
##                                  Homogenity of Variance
##                                  (With Correction Factor)
## 
## Estimated Parameter(s):          tratamiento.1 = 0.0003130000
##                                  tratamiento.2 = 0.0016003333
##                                  tratamiento.3 = 0.0003523333
##                                  tratamiento.4 = 0.0029170000
##                                  tratamiento.5 = 0.0002813333
## 
## Data:                            y
## 
## Grouping Variable:               tratamiento
## 
## Data Source:                     dat1
## 
## Sample Sizes:                    tratamiento.1 = 3
##                                  tratamiento.2 = 3
##                                  tratamiento.3 = 3
##                                  tratamiento.4 = 3
##                                  tratamiento.5 = 3
## 
## Test Statistic:                  Chisq = 4.07296
## 
## Test Statistic Parameter:        df = 4
## 
## P-value:                         0.3962219

Evaluación de la Homogeneidad de Varianzas

El valor p obtenido es 0.3962, que es mayor que 0.05, lo que indica que no se rechaza la hipótesis nula de que las varianzas de los grupos son iguales. Esto sugiere que no hay evidencia suficiente para afirmar que las varianzas son diferentes entre los grupos; por lo tanto, las varianzas pueden considerarse homogéneas.

Usando el paquete vartest

vartest::bartletts.test(y~tratamiento, dat1, alpha = 0.05)

## 
##   Bartlett's Test (alpha = 0.05) 
## ------------------------------------------------------------- 
##   data : y and tratamiento 
## 
##   statistic  : 3.959822 
##   df         : 4 
##   p.value    : 0.4114707 
## 
##   Result     : Variances are homogeneous. 
## -------------------------------------------------------------

Prueba de Levene

usando el paquete car

# paquete car 
car::leveneTest(y~tratamiento,data=dat1)

Usando el paquete vartest

vartest::levene.test(y~tratamiento, dat1,  center="median", deviation="absolute")

## 
##   Levene's Test (alpha = 0.05) 
## ------------------------------------------------------------- 
##   data : y and tratamiento 
## 
##   statistic  : 0.7889914 
##   num df     : 4 
##   denom df   : 10 
##   p.value    : 0.5580833 
## 
##   Result     : Variances are homogeneous. 
## -------------------------------------------------------------

Usando el paquete lawstat

# paquete lawstat
lawstat::levene.test( dat1$y, dat1$tratamiento, location="mean")

## 
##  Classical Levene's test based on the absolute deviations from the mean
##  ( none not applied because the location is not set to median )
## 
## data:  dat1$y
## Test Statistic = 1.6382, p-value = 0.2397

lawstat::levene.test( dat1$y, dat1$tratamiento, location="median")

## 
##  Modified robust Brown-Forsythe Levene-type test based on the absolute
##  deviations from the median
## 
## data:  dat1$y
## Test Statistic = 0.78899, p-value = 0.5581

Usando el paquete EnvStats

EnvStats::varGroupTest( y~tratamiento, data = dat1, test = "Levene")

## 
## Results of Hypothesis Test
## --------------------------
## 
## Null Hypothesis:                 Ratio of each pair of variances = 1
## 
## Alternative Hypothesis:          At least one variance differs
## 
## Test Name:                       Levene's Test for
##                                  Homogenity of Variance
## 
## Estimated Parameter(s):          tratamiento.1 = 0.0003130000
##                                  tratamiento.2 = 0.0016003333
##                                  tratamiento.3 = 0.0003523333
##                                  tratamiento.4 = 0.0029170000
##                                  tratamiento.5 = 0.0002813333
## 
## Data:                            y
## 
## Grouping Variable:               tratamiento
## 
## Data Source:                     dat1
## 
## Sample Sizes:                    tratamiento.1 = 3
##                                  tratamiento.2 = 3
##                                  tratamiento.3 = 3
##                                  tratamiento.4 = 3
##                                  tratamiento.5 = 3
## 
## Test Statistic:                  F = 1.63824
## 
## Test Statistic Parameters:       num df   =  4
##                                  denom df = 10
## 
## P-value:                         0.2397448