Paquetes necesario para la lectura de archivos excel

a) Identificación de Factores

Factores del Diseño de Cuadrado Latino (DCL)

  1. Factor de Tratamiento:
    • Ingredientes: A, B, C, D, E
  2. Factor de Fila:
    • Día: Representa los diferentes días en los que se realizan los experimentos.
  3. Factor de Columna:
    • Lote: Representa los diferentes lotes de material utilizados en el experimento.

b) Hipótesis Estadísticas y Modelo Estadístico

Hipótesis Estadísticas

  • Hipótesis Nula (H0): No hay efecto de los ingredientes sobre el tiempo de reacción (los tratamientos son iguales).
  • Hipótesis Alternativa (H1): Hay un efecto de los ingredientes sobre el tiempo de reacción (al menos un tratamiento es diferente).

Modelo Estadístico para un Diseño de Cuadrado Latino (DCL)

El modelo estadístico para un Diseño de Cuadrado Latino se expresa de la siguiente manera:

\[ Y_{ijk} = \mu + \tau_i + \rho_j + \lambda_k + \epsilon_{ijk} \]

Donde:

  • \(Y_{ijk}\): Tiempo de reacción para el tratamiento \(i\) en la fila \(j\) y columna \(k\).
  • \(\mu\): Media general del tiempo de reacción.
  • \(\tau_i\): Efecto del \(i\)-ésimo tratamiento (Ingrediente).
  • \(\rho_j\): Efecto del \(j\)-ésimo día (Fila).
  • \(\lambda_k\): Efecto del \(k\)-ésimo lote (Columna).
  • \(\epsilon_{ijk}\): Error aleatorio asociado a la observación.
library(readxl)
## Warning: package 'readxl' was built under R version 4.4.2
dat1=read_excel("ejemplo_dcl.xlsx", sheet = "dcl-dia")
head(dat1)

`

print(dat1, n = Inf) 
## # A tibble: 25 × 5
##       id  fila columna tratamiento     y
##    <dbl> <dbl>   <dbl> <chr>       <dbl>
##  1     1     1       1 A               8
##  2     2     1       2 B               7
##  3     3     1       3 D               1
##  4     4     1       4 C               7
##  5     5     1       5 E               3
##  6     6     2       1 C              11
##  7     7     2       2 E               2
##  8     8     2       3 A               7
##  9     9     2       4 D               3
## 10    10     2       5 B               8
## 11    11     3       1 B               4
## 12    12     3       2 A               9
## 13    13     3       3 C              10
## 14    14     3       4 E               1
## 15    15     3       5 D               5
## 16    16     4       1 D               6
## 17    17     4       2 C               8
## 18    18     4       3 E               6
## 19    19     4       4 B               6
## 20    20     4       5 A              10
## 21    21     5       1 E               4
## 22    22     5       2 D               2
## 23    23     5       3 B               3
## 24    24     5       4 A               8
## 25    25     5       5 C               8
print(dat1, max = Inf)
## # A tibble: 25 × 5
##       id  fila columna tratamiento     y
##    <dbl> <dbl>   <dbl> <chr>       <dbl>
##  1     1     1       1 A               8
##  2     2     1       2 B               7
##  3     3     1       3 D               1
##  4     4     1       4 C               7
##  5     5     1       5 E               3
##  6     6     2       1 C              11
##  7     7     2       2 E               2
##  8     8     2       3 A               7
##  9     9     2       4 D               3
## 10    10     2       5 B               8
## # ℹ 15 more rows
dat1$fila=as.factor(dat1$fila)
dat1$columna=as.factor(dat1$columna)
dat1$tratamiento = as.factor(dat1$tratamiento)

ESTADISTICA DESCRIPTIVA

Grafico de box-plot

boxplot(y ~ tratamiento, data = dat1, xlab = "tratamiento",  ylab = "reaccion proceso-quimico", las=1, 
        cex.axis=1.2, cex.lab=1.2)

DCL

ANOVA

fit1 = lm(y ~ tratamiento + fila + columna, data=dat1)
anova(fit1)

Efecto del Tratamiento

El tratamiento tiene un valor p de 0.0004877, que es significativamente menor que (\(\alpha = 0.05\)). Esto indica que hay diferencias significativas en la variable de respuesta \(y\) entre al menos uno de los tratamientos. Por lo tanto, rechazamos la hipótesis nula que establece que todos los tratamientos tienen el mismo efecto sobre \(y\). Esto sugiere que al menos uno de los métodos de tratamiento es efectivo para reducir el conteo de partículas.

Efecto de la Fila

El efecto de la fila (día) tiene un valor p de 0.3476182, que es mucho mayor que (\(\alpha = 0.05\)). Esto indica que no hay evidencia suficiente para afirmar que las diferencias entre los días tienen un impacto significativo en la variable de respuesta \(y\). Por lo tanto, no rechazamos la hipótesis nula para el efecto de la fila.

Efecto de la Columna

El efecto de la columna (lote) también muestra un valor p de 0.4550143, que es igualmente mayor que (\(\alpha = 0.05\)). Esto sugiere que no hay diferencias significativas entre los lotes en relación con la variable de respuesta \(y\). Por lo tanto, no rechazamos la hipótesis nula para el efecto de la columna.

Conclusiones Generales

Los grados de libertad (gl) para cada uno delas sumas de los cuadrados (SC) están en el vector anova(fit1)$Df

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

Los gl para los tratamientos, los factores fila, columna y el error son respectivamente

# tratamientos
gl.trat=gl_anova[1]; gl.trat
## [1] 4
# factor fila
gl.fil=gl_anova[2]; gl.fil
## [1] 4
# factor columna
gl.col=gl_anova[3]; gl.col
## [1] 4
# error
gl.error=gl_anova[4]; gl.error
## [1] 12

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

sc_anova=anova(fit1)$Sum

Luego, la SC para los tratamientos, los factores fila, columna y el error son respectivamente

sctr=sc_anova[1];sctr 
## [1] 141.44
scfil=sc_anova[2];scfil
## [1] 15.44
sccol=sc_anova[3];sccol
## [1] 12.24
sce=sc_anova[4];sce
## [1] 37.52
sct=sum(sc_anova);sct
## [1] 206.64

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

cm_anova=anova(fit1)$Mean; cm_anova
## [1] 35.360000  3.860000  3.060000  3.126667

Luego, los CM para los tratamientos, los factores fila, columna y el error son respectivamente

cmtrat=cm_anova[1]; cmtrat
## [1] 35.36
cmfil=cm_anova[2]; cmfil
## [1] 3.86
cmcol=cm_anova[3]; cmcol
## [1] 3.06
cme=cm_anova[4]; cme
## [1] 3.126667

d) 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_dcl1=LSD.test(dat1$y, dat1$tratamiento, DFerror=gl.error, MSerror=cme, p.adj = "none")
lsd_dcl1
## $statistics
##    MSerror Df Mean       CV  t.value      LSD
##   3.126667 12 5.88 30.07208 2.178813 2.436636
## 
## $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 Q50 Q75
## A    8.4 1.140175 5 0.7907802 6.677038 10.122962   7  10   8   8   9
## B    5.6 2.073644 5 0.7907802 3.877038  7.322962   3   8   4   6   7
## C    8.8 1.643168 5 0.7907802 7.077038 10.522962   7  11   8   8  10
## D    3.4 2.073644 5 0.7907802 1.677038  5.122962   1   6   2   3   5
## E    3.2 1.923538 5 0.7907802 1.477038  4.922962   1   6   2   3   4
## 
## $comparison
## NULL
## 
## $groups
##   dat1$y groups
## C    8.8      a
## A    8.4      a
## B    5.6      b
## D    3.4      b
## E    3.2      b
## 
## attr(,"class")
## [1] "group"
plot(lsd_dcl1, ,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
## B-A -2.8 -6.3646078  0.7646078 0.1539433
## C-A  0.4 -3.1646078  3.9646078 0.9960012
## D-A -5.0 -8.5646078 -1.4353922 0.0055862
## E-A -5.2 -8.7646078 -1.6353922 0.0041431
## C-B  3.2 -0.3646078  6.7646078 0.0864353
## D-B -2.2 -5.7646078  1.3646078 0.3365811
## E-B -2.4 -5.9646078  1.1646078 0.2631551
## D-C -5.4 -8.9646078 -1.8353922 0.0030822
## E-C -5.6 -9.1646078 -2.0353922 0.0023007
## E-D -0.2 -3.7646078  3.3646078 0.9997349

Grafico para las comparaciones múltiples

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

Interpretación de Comparaciones Múltiples

Fisher-LSD

  • LSD: 2.436636 (diferencias significativas si exceden este valor).
  • Medias de Tratamientos:
    • A: 8.4
    • B: 5.6
    • C: 8.8
    • D: 3.4
    • E: 3.2
  • Grupos:
    • Grupo A: C y A son similares.
    • Grupo B: B, D y E son significativamente diferentes de A y

Conclusión: A y C son efectivos; B, D y E son menos efectivos.

Tukey

  • Diferencias Significativas:
    • D - A: -5.0 (p = 0.0055862)
    • E - A: -5.2 (p = 0.0041431)
    • D - C: -5.4 (p = 0.0030822)
    • E - C: -5.6 (p = 0.0023007)
  • No Significativas: Comparaciones entre B y A/C.

Conclusión: D y E son inferiores a A y C; B no muestra diferencias significativas.

Conclusión General

Los análisis realizados indican que los tratamientos A y C son significativamente más efectivos en la reducción del conteo de partículas en comparación con los tratamientos B, D y E. Esto sugiere que para lograr una mayor eficacia en el proceso de reducción de partículas, se deben priorizar los métodos A y C, mientras que los métodos B, D y E requieren reconsideración o modificación para mejorar su rendimiento.

VERIFICACION DE HIPOTESIS

e) Análisis de 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   8      8.44     -0.44
## 2   7      4.64      2.36
## 3   1      2.24     -1.24
## 4   7      7.24     -0.24
## 5   3      3.44     -0.44
## 6  11      9.84      1.16
## 7   2      3.24     -1.24
## 8   7      8.24     -1.24
## 9   3      2.84      0.16
## 10  8      6.84      1.16
## 11  4      6.24     -2.24
## 12  9      8.04      0.96
## 13 10      8.24      1.76
## 14  1      2.24     -1.24
## 15  5      4.24      0.76
## 16  6      5.44      0.56
## 17  8      9.84     -1.84
## 18  6      4.04      1.96
## 19  6      6.04     -0.04
## 20 10     10.64     -0.64
## 21  4      3.04      0.96
## 22  2      2.24     -0.24
## 23  3      4.24     -1.24
## 24  8      6.64      1.36
## 25  8      8.84     -0.84

Analisis

El análisis de residuos revela que algunos puntos de datos presentan residuos significativos, lo que indica que el modelo no está capturando adecuadamente todas las variaciones en los datos. En particular, se observan subestimaciones y sobreestimaciones en ciertos casos, lo que sugiere que el modelo podría no estar considerando todas las variables relevantes o que la relación entre las variables no es lineal.

conclusion

En general, el modelo presenta un ajuste razonable, pero es necesario realizar ajustes adicionales o considerar la inclusión de variables adicionales para mejorar su precisión. Un análisis más profundo de los residuos es fundamental para verificar el cumplimiento de los supuestos del modelo y asegurar un ajuste adecuado a los datos.

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
Se puede observar que ningun punto sale de la banda de confianza, esto indica que los errores son normales

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

Considerando como nivel de significancia \(\alpha\)=0.05, usando la prueba SW se tiene que el valor p=0.5476 que es mayor que 0.05, asi los errores son normales.

shapiro.test(fit1$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  fit1$residuals
## W = 0.96606, p-value = 0.5476
ks.test(fit1$residuals,"pnorm")
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  fit1$residuals
## D = 0.17251, p-value = 0.4008
## alternative hypothesis: two-sided

Independencia de los errores

Se observa que no hay patron aparente en los residuos. Po esta razón se puede indicar que existe la independencia de los errores

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

Homogeneidad de varianzas (Homocedasticidad)

Del grafico de los residuos vs los valores ajustados se puede observar que no hay un patron aparente, así se puede indicar considerando el grafico que la homogeneidad de las varianzas se cumple

#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 = 1.5544, df = 4, p-value = 0.817

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):          A = 1.3
##                                  B = 4.3
##                                  C = 2.7
##                                  D = 4.3
##                                  E = 3.7
## 
## Data:                            y
## 
## Grouping Variable:               tratamiento
## 
## Data Source:                     dat1
## 
## Sample Sizes:                    A = 5
##                                  B = 5
##                                  C = 5
##                                  D = 5
##                                  E = 5
## 
## Test Statistic:                  Chisq = 1.578303
## 
## Test Statistic Parameter:        df = 4
## 
## P-value:                         0.8126864

Usando el paquete vartest

vartest::bartletts.test(y~tratamiento, dat1, alpha = 0.05)
## 
##   Bartlett's Test (alpha = 0.05) 
## ------------------------------------------------------------- 
##   data : y and tratamiento 
## 
##   statistic  : 1.554389 
##   df         : 4 
##   p.value    : 0.8169654 
## 
##   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.4444444 
##   num df     : 4 
##   denom df   : 20 
##   p.value    : 0.7751009 
## 
##   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 = 0.79715, p-value = 0.541
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.44444, p-value = 0.7751

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):          A = 1.3
##                                  B = 4.3
##                                  C = 2.7
##                                  D = 4.3
##                                  E = 3.7
## 
## Data:                            y
## 
## Grouping Variable:               tratamiento
## 
## Data Source:                     dat1
## 
## Sample Sizes:                    A = 5
##                                  B = 5
##                                  C = 5
##                                  D = 5
##                                  E = 5
## 
## Test Statistic:                  F = 0.797153
## 
## Test Statistic Parameters:       num df   =  4
##                                  denom df = 20
## 
## P-value:                         0.5410408