Laboratorio #2 - Implementacion Diseño Factorial

Asignacion de los posibles valores que pueden tomar las variables

library(dplyr)
valores_ABC <- c(-1,-1,-1,
                 1,-1,-1,
                -1,1,-1,
                 1,1,-1,
                -1,-1,1,
                 1,-1,1,
                -1,1,1,
                 1,1,1,
                -1,-1,-1,
                 1,-1,-1,
                -1,1,-1,
                 1,1,-1,
                -1,-1,1,
                 1,-1,1,
                -1,1,1,
                 1,1,1)
valores_ABC <- matrix(ABC_values, ncol = 3, byrow = TRUE)
valores_ABC

      [,1] [,2] [,3]
 [1,]   -1    1   -1
 [2,]    1   -1    1
 [3,]   -1    1   -1
 [4,]    1   -1    1
 [5,]   -1    1   -1
 [6,]    1   -1   -1
 [7,]    1    1   -1
 [8,]   -1    1    1
 [9,]   -1   -1    1
[10,]    1   -1   -1
[11,]    1    1   -1
[12,]   -1   -1   -1
[13,]    1    1    1
[14,]    1   -1   -1
[15,]   -1   -1    1
[16,]    1    1    1

Asignacion de valores a las replicas del experimento

rep_exp <- c(247,470,429,435,837,551,775,660,400,446,405,445,850,670,865,530)
rep_exp

 [1] 247 470 429 435 837 551 775 660 400 446 405 445 850 670 865 530

plasma <- data.frame(cbind(A = valores_ABC[,1], B = valores_ABC[,2], C = valores_ABC[,3],  replica = rep_exp))
head(plasma)

Valores totales asignados a las replicas

val_total <- c(647,916,834,880,1687,1221,1640,1190)
val_total <- matrix(val_total)
rownames(val_total) <- c("(1)","a","b","ab","c","ac","bc","abc")
val_total

    [,1]
(1)  647
a    916
b    834
ab   880
c   1687
ac  1221
bc  1640
abc 1190

El efecto promedio de A, donde \(n\) es la cantidad de replicas. para este laboratorio son 2. Las cantidades dentro de los corchetes \([\hspace{0.1cm}]\) son los contrastes. \[A = \frac{1}{4n}[a + ab + ac + abc - (1) - b - c - bc]\]

\[B = \frac{1}{4n}[b + ab + bc + abc - (1) - a - c - ac]\]

\[C = \frac{1}{4n}[c + ac + bc + abc - (1) - a - c - ab]\]

Calculamos los efectos #### A

numero_rep <- 2
A_cont <- val_total["a",] + val_total["ab",] + val_total["ac",] + val_total["abc",] - val_total["(1)",] - val_total["b",] - val_total["c",] - val_total["bc",]
A_cont

   a 
-601 

A_effecto <- (1/(4*rep_num)) * A_cont
A_effecto

      a 
-75.125 

B

B_cont <- val_total["b",] + val_total["ab",] + val_total["bc",] + val_total["abc",] - val_total["(1)",] - val_total["a",] - val_total["c",] - val_total["ac",]
B_cont

 b 
73 

B_effecto <- (1/(4*rep_num)) * B_cont
B_effecto

    b 
9.125 

C

C_cont <- val_total["c",] + val_total["ac",] + val_total["bc",] + val_total["abc",] - val_total["(1)",] - val_total["a",] - val_total["b",] - val_total["ab",]
C_cont

   c 
2461 

C_effecto <- (1/(4*rep_num)) * C_cont
C_effecto

      c 
307.625 

Interacion AB es la diferencia en los promedios entre las corridas de dos planos diagonales \[AB = \frac{1}{4n}[abc + ab + c + (1) - bc - b - ac - a]\]

\[AC = \frac{1}{4n}[(1) + b + ac + abc - a - ab - c - bc]\]

\[BC = \frac{1}{4n}[(1) + a + bc + abc - b - ab - c - ac]\]

AB_cont <- val_total["abc",] + val_total["ab",] + val_total["c",] + val_total["(1)",] - val_total["bc",] - val_total["b",] - val_total["ac",] - val_total["a",]
AB_cont

 abc 
-207 

AB_effecto <- (1/(4*rep_num)) * AB_cont
AB_effecto

    abc 
-25.875 

AC_cont <- val_total["(1)",] + val_total["b",] + val_total["ac",] + val_total["abc",] - val_total["a",] - val_total["ab",] - val_total["c",] - val_total["bc",]
AC_cont

  (1) 
-1231 

AC_effecto <- (1/(4*rep_num)) * AC_cont
AC_effecto

     (1) 
-153.875 

BC_cont <- val_total["(1)",] + val_total["a",] + val_total["bc",] + val_total["abc",] - val_total["b",] - val_total["ab",] - val_total["c",] - val_total["ac",]
BC_cont

 (1) 
-229 

BC_effecto <- (1/(4*rep_num)) * BC_cont
BC_effecto

    (1) 
-28.625 

Interaccion ABC es la diferencia de dos promedios, AB - C

ABC_cont <- val_total["abc",] + val_total["c",] + val_total["b",] + val_total["a",] - val_total["bc",] - val_total["ac",] - val_total["ab",] - val_total["(1)",]
ABC_cont

abc 
239 

ABC_effecto <- (1/(4*rep_num)) * ABC_cont
ABC_effecto

   abc 
29.875 

calculamos la suma de los cuadrados de los efectos de la siguiente manera \[SS = \frac{(contrast)^2}{8n}\]

A_ss <- (A_cont^2)/(8*numero_rep)
B_ss <- (B_cont^2)/(8*numero_rep)
C_ss <- (C_cont^2)/(8*numero_rep)
AB_ss <- (AB_cont^2)/(8*numero_rep)
AC_ss <- (AC_cont^2)/(8*numero_rep)
BC_ss <- (BC_cont^2)/(8*numero_rep)
ABC_ss <- (ABC_cont^2)/(8*numero_rep)
c(A_ss,B_ss,C_ss,AB_ss,AC_ss,BC_ss,ABC_ss)

          a           b           c         abc         (1)         (1)         abc 
 22575.0625    333.0625 378532.5625   2678.0625  94710.0625   3277.5625   3570.0625 

Analizaremos la varianza, encontraremos los P-values y se encontraremos los factores estadisticamente significativos. Utilizaremos anova() Analysis of Variance, la cual recibe una regresion, y por ende utilizaremos todos los efectos: A, B, C, AB, AC, BC y ABC.

fit_AOV <- lm(replica ~ A + B + C + (A*B) + (A*C) + (B*C) + (A*B*C), data = plasma)
anova(fit_AOV)

Analysis of Variance Table

Response: replica
          Df Sum Sq Mean Sq F value Pr(>F)
A          1    947     947  0.0174 0.8984
B          1  13821   13821  0.2536 0.6281
C          1  22325   22325  0.4096 0.5401
A:B        1  13271   13271  0.2435 0.6350
A:C        1  28377   28377  0.5206 0.4911
B:C        1  11038   11038  0.2025 0.6646
A:B:C      1  11834   11834  0.2171 0.6537
Residuals  8 436058   54507               

En Sum Sq se observa los mismos datos obtenidos en la suma de cuadrados. El error es de \(31,996\). Las variables mas significativas para la regresion son: A, C y AC.

plasma_fit <- lm(replica ~ A + C + (A*C), data = plasma)
summary(plasma_fit)

Call:
lm(formula = replica ~ A + C + (A * C), data = plasma)

Residuals:
    Min      1Q  Median      3Q     Max 
-242.50 -126.61  -42.88  126.85  347.50 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  567.954     51.720  10.981 1.29e-07 ***
A              2.371     51.720   0.046    0.964    
C             38.504     51.720   0.744    0.471    
A:C          -37.579     51.720  -0.727    0.481    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 203.5 on 12 degrees of freedom
Multiple R-squared:  0.0756,    Adjusted R-squared:  -0.1555 
F-statistic: 0.3271 on 3 and 12 DF,  p-value: 0.8058

plasma_fit$coefficients

(Intercept)           A           C         A:C 
 567.954167    2.370833   38.504167  -37.579167 

Otra forma para determinar las variables mas significativas es calculando el error estandar de los efectos \[Var(Effect) = Var\bigg(\frac{Contrast}{n2^{k-1}}\bigg) = \frac{1}{(n2^{k-1})^2}Var(Contrast)\]

Cada constraste es una combinacion lineal de los 2K totales de los tratamientos, y cada total consta de n observaciones \[Var(Contrast) = n2^k\sigma^2\]

se puede reescribir la ecuacion del error como \[Var(Effect)=\frac{1}{(n2^{k-1})^2}n2^k\sigma^2 = \frac{1}{n2^{k-2}}\sigma^2\]

El error estimado se encuentra como: \[se(Effect) =\sqrt{\frac{1}{n2^{k-2}}s^2}\]

el valor estimado del error de minimos cuadrados se obtiene de la columna mean sq, por lo que:

MSE <- anova(fit_AOV)$`Mean Sq`[8]
MSE

[1] 54507.29

se_effecto <- sqrt(1/(numero_rep*(2^(3-2)))*(MSE))
se_effecto

[1] 116.734

para calcular los limites de los errores en el estimado de los efectos con dos desviaciones estandar son:

lim <- c(A_effecto - (2*se_effecto), A_effecto + (2*se_effecto),
            B_effecto - (2*se_effecto), B_effecto + (2*se_effecto),
            C_effecto - (2*se_effecto), C_effecto + (2*se_effecto),
            AB_effecto - (2*se_effecto), AB_effecto + (2*se_effecto),
            AC_effecto - (2*se_effecto), AC_effecto + (2*se_effecto),
            BC_effecto - (2*se_effecto), BC_effecto + (2*se_effecto),
            ABC_effecto - (2*se_effecto), ABC_effecto + (2*se_effecto))
lim <- matrix(lim, ncol = 2, byrow = TRUE)
rownames(lim) <- c("A","B","C","AB","AC","BC","ABC")
colnames(lim) <- c("low", "high")
lim

           low      high
A   -308.59297 158.34297
B   -224.34297 242.59297
C     74.15703 541.09297
AB  -259.34297 207.59297
AC  -387.34297  79.59297
BC  -262.09297 204.84297
ABC -203.59297 263.34297