R Notebook

Ejemplo de un experimento \(3^2\)

(Planear) el experimento

La planeación de las pruebas experimentales (corridas) requiere comprender el objeto o sistema bajo estudio, tal que permita plantear un sistema de preguntas de investigación e hipótesis, al haber identificado con una base teórica o práctica las variables candidato a factores de control (posibles causas especiales de variación), así como los rangos de operación en que los ajustes de los factores producen la respuesta deseada. Este conocimiento se resume en la lista de variables y sus niveles de control para el experimento.

library (DoE.base)
myvar1 <- list (ra=c(140,160,180),
                ce=c(200,250,300))
t0 <- fac.design(factor.names = myvar1, randomize = F, replications = 2)

creating full factorial with 9 runs ...

t0

class=design, type= full factorial 
NOTE: columns run.no and run.no.std.rp  are annotation,  not part of the data frame

(Experimentar) y recolectar los datos

La realización del experimento require la ejecución cuidadosa y controlada de las corridas experimentales con base en el plan experimental. En esta etapa los errores de medición (incertidumbre, exactitud y precisión) podrían incrementar el error hasta impedir que las variaciones de los factores de control se puedan observar (\(\beta\) probabilidad de un falso positivo), así como también errores que provoquen sesgos.

Y <- c(60.75, 184.21, 296.09,  87.64, 267.71, 414.72,  96.29, 333.75, 521.53,  49.72, 185.33, 292.21,  84.81, 260.80, 420.75, 112.04, 346.84, 532.73)

t0 <- add.response(t0,Y)
t0

class=design, type= full factorial 
NOTE: columns run.no and run.no.std.rp  are annotation,  not part of the data frame

(Verificar) la hipótesis mediante el análisis estadístico y su interpretación

El modelado estadístico permitirá relacionar los resultados experimentales con la hipótesis original y verificarla.

Nota: En este caso, los coeficientes de la función {DoE.base::halfnormal()} no corresponden a los coeficientes de la regresión, por que se trata de un experimento con tres niveles.

plot (t0) #Resumen gráfico: ra tiene el efecto mayor

The first four factors were selected. Use argument select for choosing what is plotted.

hf <- halfnormal(t0, alpha = .05, method = "LW98") #Ra es el factor significativo

Significant effects (alpha=0.05, LW98 method):
[1] ra.L      ce.L      ra.L:ce.L ra.Q      

hf$signif #Factores significativos

[1] "ra.Q"      "ra.L:ce.L" "ce.L"      "ra.L"     

hf$coef

        ra.L         ra.Q         ce.L         ce.Q    ra.L:ce.L    ra.Q:ce.L    ra.L:ce.Q    ra.Q:ce.Q 
135.18325642  -7.38533749  59.52736365  -2.41084129  30.67500000  -2.80399781  -0.10969655  -1.13055556 
          e1           e2           e3           e4           e5           e6           e7           e8 
 -0.03820365  -2.80995011   1.40744393  -0.89460327   1.22018323  -0.64739953   0.18666667  -1.83833333 
          e9 
  2.35059171 

hf$LCs

named list()

fit <- lm (Y~ra+ce+ra:ce, t0)

anova (fit)

Analysis of Variance Table

Response: Y
          Df Sum Sq Mean Sq F value    Pr(>F)    
ra         2 329923  164962 3831.67 6.485e-14 ***
ce         2  63888   31944  741.98 1.025e-10 ***
ra:ce      4  17102    4275   99.31 1.936e-07 ***
Residuals  9    387      43                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

summary (fit)

Call:
lm.default(formula = Y ~ ra + ce + ra:ce, data = t0)

Residuals:
   Min     1Q Median     3Q    Max 
-7.875 -3.345  0.000  3.345  7.875 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 252.6622     1.5465 163.373  < 2e-16 ***
ra.L        234.1443     2.6787  87.410 1.70e-14 ***
ra.Q        -12.7918     2.6787  -4.775  0.00101 ** 
ce.L        103.1044     2.6787  38.491 2.68e-11 ***
ce.Q         -4.1757     2.6787  -1.559  0.15346    
ra.L:ce.L    92.0250     4.6396  19.835 9.77e-09 ***
ra.Q:ce.L    -8.4120     4.6396  -1.813  0.10323    
ra.L:ce.Q    -0.3291     4.6396  -0.071  0.94500    
ra.Q:ce.Q    -3.3917     4.6396  -0.731  0.48336    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.561 on 9 degrees of freedom
Multiple R-squared:  0.9991,    Adjusted R-squared:  0.9982 
F-statistic:  1193 on 8 and 9 DF,  p-value: 1.078e-12

library (sjPlot)
plot_model(fit, type = "est", show.values = T) #Forest plot efectos

plot_model(fit, type = "int") #Interacciones

plot_model(fit, type = "diag") #Diagnóstico del modelo

[[1]]

[[2]]

`geom_smooth()` using formula 'y ~ x'

[[3]]

[[4]]

`geom_smooth()` using formula 'y ~ x'

(Aplicar) Trasladar los resultados del análisis al problema original.

El resultado que se traslada a la solución del problema es la tabla de los intervalos de confianza con los efectos y la gráfica de efectos de la interacción, que es una ayuda para interpretar los valores estimados de distancia con las combinaciones de RA y CE.

library (effects)
ef <- allEffects(fit)
summary (ef)

 model: Y ~ ra + ce + ra:ce

 ra*ce effect
     ce
ra        200     250     300
  140  55.235  86.225 104.165
  160 184.770 264.255 340.295
  180 294.150 417.735 527.130

 Lower 95 Percent Confidence Limits
     ce
ra          200       250       300
  140  44.73945  75.72945  93.66945
  160 174.27445 253.75945 329.79945
  180 283.65445 407.23945 516.63445

 Upper 95 Percent Confidence Limits
     ce
ra          200       250      300
  140  65.73055  96.72055 114.6605
  160 195.26555 274.75055 350.7905
  180 304.64555 428.23055 537.6255

as.data.frame(ef)

$`ra:ce`

plot (ef)

NA
NA

Convenience resources in {DoE.base}

There is DoE.base::lm function that can get an automatic statistical model directly from the experimental table with a response column; but, caution should be excercise because the generated formula is not necessarily the best for statistics analysis. In this case there are four non significant factors included.

lm1 <- lm(t0)
summary(lm1)

Number of observations used: 18 
Formula:
Y ~ (ra + ce)^2

Call:
lm.default(formula = fo, data = model.frame(fo, data = formula))

Residuals:
   Min     1Q Median     3Q    Max 
-7.875 -3.345  0.000  3.345  7.875 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 252.6622     1.5465 163.373  < 2e-16 ***
ra.L        234.1443     2.6787  87.410 1.70e-14 ***
ra.Q        -12.7918     2.6787  -4.775  0.00101 ** 
ce.L        103.1044     2.6787  38.491 2.68e-11 ***
ce.Q         -4.1757     2.6787  -1.559  0.15346    
ra.L:ce.L    92.0250     4.6396  19.835 9.77e-09 ***
ra.Q:ce.L    -8.4120     4.6396  -1.813  0.10323    
ra.L:ce.Q    -0.3291     4.6396  -0.071  0.94500    
ra.Q:ce.Q    -3.3917     4.6396  -0.731  0.48336    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.561 on 9 degrees of freedom
Multiple R-squared:  0.9991,    Adjusted R-squared:  0.9982 
F-statistic:  1193 on 8 and 9 DF,  p-value: 1.078e-12