R Notebook

Ejemplo de un experimento \(2^5\)

(PLAN) 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.

# EXPERIMENTO 2^5 ----
# https://sigmazone.com/catapult-grid/

library (DoE.base)

Warning: package ‘DoE.base’ was built under R version 4.1.3
Loading required package: grid
Loading required package: conf.design

Attaching package: ‘DoE.base’

The following objects are masked from ‘package:stats’:

    aov, lm

The following object is masked from ‘package:graphics’:

    plot.design

The following object is masked from ‘package:base’:

    lengths

vars <- list(ra=c(140,180),
             fa=c(120,135),
             ce=c(200,300),
             pe=c(150,200),
             bp=c(150,250)
)
tabla <- fac.design(factor.names = vars, randomize = F)

creating full factorial with 32 runs ...

tabla

class=design, type= full factorial 

(DO) Realizar el experimento 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(76.68,214.35,0,155.03,142.48,378.26,15.68,290.77,
       105.13,276.46,5.65,211.87,184.94,493.15,37.23,
       378.04,102.74,257.81,7.01,190.75,184.82,462.52,
       26.59,360.94,131.55,333.75,15.22,255.12,226.40,
       602.26,50.29,465.19)

tab.1 <- add.response(tabla,Y)

tab.1

class=design, type= full factorial 

(CHECK) Análisis estadístico y su interpretación

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

plot (tab.1) #Resumen gráfico: ra tiene el efecto mayor

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

hf <- halfnormal(tab.1, alpha = .05) #Ra es el factor significativo

Significant effects (alpha=0.05, Lenth method):
 [1] ra1     ce1     fa1     ra1:ce1 pe1     bp1     ra1:pe1 fa1:ce1 ra1:bp1 ra1:fa1 ce1:pe1 ce1:bp1 

hf$signif #Factores significativos

 [1] "ce1:bp1" "ce1:pe1" "ra1:fa1" "ra1:bp1" "fa1:ce1" "ra1:pe1" "bp1"     "pe1"     "ra1:ce1" "fa1"    
[11] "ce1"     "ra1"    

hf$coef

        ra1         fa1         ce1         pe1         bp1     ra1:fa1     ra1:ce1     ra1:pe1     ra1:bp1 
125.4331250 -53.3725000  61.2637500  28.3068750  22.1012500   8.9443750  34.7356250  15.7812500  11.0493750 
    fa1:ce1     fa1:pe1     fa1:bp1     ce1:pe1     ce1:bp1     pe1:bp1        lof1        lof2        lof3 
-12.2587500  -5.0668750  -4.7987500   7.6581250   6.5525000   2.1056250  -1.4822612   0.7218978   0.3922572 
       lof4        lof5        lof6        lof7        lof8        lof9       lof10       lof11       lof12 
  2.5498200  -0.3259142  -0.6462058   2.9371875  -0.9691526  -2.7894393  -1.0513449  -2.3464270  -0.7771393 
      lof13       lof14       lof15       lof16 
 -1.7691652  -1.7830390   2.2324436   1.2942026 

FrF2::DanielPlot(tab.1)

fit <- lm (Y~ra+ce+fa+ra:ce+pe+bp+ra:pe+ce:fa+ra:bp+ra:fa, tab.1)

anova (fit)

Analysis of Variance Table

Response: Y
          Df Sum Sq Mean Sq   F value    Pr(>F)    
ra         1 503471  503471 1634.4866 < 2.2e-16 ***
ce         1 120104  120104  389.9097 4.842e-15 ***
fa         1  91156   91156  295.9320 7.454e-14 ***
pe         1  25641   25641   83.2417 9.413e-09 ***
bp         1  15631   15631   50.7447 5.022e-07 ***
ra:ce      1  38610   38610  125.3450 2.588e-10 ***
ra:pe      1   7970    7970   25.8726 4.888e-05 ***
ce:fa      1   4809    4809   15.6117 0.0007301 ***
ra:bp      1   3907    3907   12.6833 0.0018447 ** 
ra:fa      1   2560    2560    8.3111 0.0089028 ** 
Residuals 21   6469     308                        
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

summary (fit)

Call:
lm.default(formula = Y ~ ra + ce + fa + ra:ce + pe + bp + ra:pe + 
    ce:fa + ra:bp + ra:fa, data = tab.1)

Residuals:
    Min      1Q  Median      3Q     Max 
-30.079  -5.452  -1.621   3.074  39.443 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  207.459      3.103  66.867  < 2e-16 ***
ra1          125.433      3.103  40.429  < 2e-16 ***
ce1           61.264      3.103  19.746 4.84e-15 ***
fa1          -53.372      3.103 -17.203 7.45e-14 ***
pe1           28.307      3.103   9.124 9.41e-09 ***
bp1           22.101      3.103   7.124 5.02e-07 ***
ra1:ce1       34.736      3.103  11.196 2.59e-10 ***
ra1:pe1       15.781      3.103   5.087 4.89e-05 ***
ce1:fa1      -12.259      3.103  -3.951  0.00073 ***
ra1:bp1       11.049      3.103   3.561  0.00184 ** 
ra1:fa1        8.944      3.103   2.883  0.00890 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 17.55 on 21 degrees of freedom
Multiple R-squared:  0.9921,    Adjusted R-squared:  0.9884 
F-statistic: 264.2 on 10 and 21 DF,  p-value: < 2.2e-16

library (sjPlot)

Warning: package ‘sjPlot’ was built under R version 4.1.3
Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!

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

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

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

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

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     

[[1]]

[[2]]

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

[[3]]

[[4]]

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

# EXPERIMENTO 3^K ----
# PAPER HELICOPTER
# https://www.jpl.nasa.gov/edu/learn/project/make-a-paper-mars-helicopter/

variables <- list(ra=c(140,160,180),
                  ce=c(200,250,300))
k3 <- fac.design(factor.names = variables, randomize = F, replications = 3)

creating full factorial with 9 runs ...

k3

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